Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
101,582 characters · 26 sections · 40 citation commands
Augmented Factor Models with Applications to Validating Market Risk Factors and Forecasting Bond Risk Premia
\def\spacingset#1{ {#1}} \spacingset{1}
\if11 \fi
\if01 {
} \fi
{\it Keywords:} Heavy tails, Forecasts; Principal components; identification.
\spacingset{1.45}
In this paper, we study the identification and estimations of factor models augmented by a set of additional covariates that are common to all individuals. Consider the following factor model:
Here $\bfm y_t=(y_{1t},..., y_{Nt})'$ is the multivariate outcome for the $t^{th}$ observation in the sample; $\bfm f_t$ is the $K$-dimensional vector of latent factors; $\ensuremath{\boldsymbol{\Lambda}}=(\ensuremath{\boldsymbol{\lambda}}_1,....,\ensuremath{\boldsymbol{\lambda}}_N)'$ is an $N\times K$ matrix of nonrandom factor loadings; $\bfm u_t=(u_{1t},...,u_{Nt})'$ denotes the vector of idiosyncratic errors. In addition to $\{\bfm y_t\}_{t=1}^T$, we also observe variables, denoted by $\bfm x_t$, that have some explanatory power on the unknown factors and hence impact on observed vector $\bfm y_t$. We model $\bfm f_t$ by using the model
for some (nonparametric) function $\bfm g = E(\bfm f_t|\bfm x_t)$. Here $\bfm g(\bfm x_t)$ is interpreted as the component of the factors that can be explained by the covariates, and $\bfsym \gamma_t$ is the components that cannot be explained by the covariates. We aim to provide an improved estimation procedure when the factors can be partially explained by several observed variables $\bfm x_t$. In addition, by accurately estimating $\bfsym \gamma_t$, we can estimate the percentage of both explained and unexplained components in the factors, which describes the proxy/explanatory power of covariates.
Note that model ((ref)) implies:
where $\operatorname{cov}(\bfm y_t)$ and $\operatorname{cov}(\bfm u_t)$ respectively denote the $N\times N$ variance-covariance matrices of $\bfm y_t$ and $\bfm u_t$; $\operatorname{cov}(\bfm f_t)$ denotes the $K\times K$ variance-covariance matrix of $\bfm f_t$. Under usual factor models without covariates, $\frac{1}{\sqrt{N}}\ensuremath{\boldsymbol{\Lambda}}$ is identified asymptotically as the first $K$ eigenvectors of $\operatorname{cov}(\bfm y_t)$ as $N\to\infty$ and can be estimated using the first $K$ eigenvectors of the sample covariance matrix of $\bfm y_t$ (e.g,, SW02, bai03).
With additional covariates, on the other hand, exact identification can be achieved through covariance of the “smoothed data”. By (ref), assuming exogeneity of $\bfm x_t$, we have $ E(\bfm y_t|\bfm x_t) = \ensuremath{\boldsymbol{\Lambda}} E(\bfm f_t | \bfm x_t) $ so that it becomes a “noiseless" factor model with smoothed data $ E(\bfm y_t|\bfm x_t)$ as input and $E(\bfm f_t | \bfm x_t)$ as latent factors. The factor loadings and latent factors can be extracted from
It is easy to see from the model that
where $\bfsym \Sigma_{f|x} = E\{E(\bfm f_t|\bfm x_t)E(\bfm f_t|\bfm x_t)'\}$ is a $K\times K$ low-dimensional positive definite matrix. This decomposition is to be compared with ((ref)), where the noise covariance $\operatorname{cov}(\bfm u_t)$ removed. Therefore, as long as $\bfsym \Sigma_{f|x}$ is of full rank, $\ensuremath{\boldsymbol{\Lambda}}$ falls in the eigenspace generated by $\bfsym \Sigma_{y|x}$. In other words, $\ensuremath{\boldsymbol{\Lambda}}$ is identifiable up to an orthogonal transformation. Because of such exact identification, we allow $N$ to be finite as a special case. The number of factors is assumed to be known throughout the paper. In practice, $K$ can be consistently estimated by many methods such as AIC, BIC-based criteria, or eigenvalue-ratio methods studied in LamYao,AH.
The above discussion prompts us the following new method to estimate the factor loadings $\ensuremath{\boldsymbol{\Lambda}}$ that incorporates the explanatory power of $\bfm x_t$: (See Section 3 for details of estimators)
(i) (robustly) regress $\{\bfm y_t\}$ on $\{\bfm x_t\}$ and obtain fitted value $\{\widehat\bfm y_t\}$;
(ii) conduct the principal components analysis (PCA) on the fitted data $(\widehat\bfm y_1,...,\widehat\bfm y_T)$ to estimate the factor loadings. \\ We employ a regression based on huber1964robust's robust M-estimation in step (i). The procedure involves a diverging truncation parameter, called adaptive Huber loss, to reduce the bias when the error distribution is asymmetric fan2017estimation. This allows our procedure to be applicable to data with heavy tails.\footnote{In this paper, by “heavy-tail" we mean tail distributions of $(\bfm u_t, \bfm y_t)$ that are heavier than the usual requirements on the high-dimensional factor model (which are either exponentially-tailed or have eighth or higher moments). But we do not allow large outliers on the covariates. }
There are two important quantities that determine the rates of convergence for the estimators: the “signal" $\bfsym \Sigma_{f|x}= E\{E(\bfm f_t|\bfm x_t)E(\bfm f_t|\bfm x_t)'\}$ and the “noise" $\operatorname{cov}(\bfsym \gamma_t)$. The rates of convergence are presented using these two quantities. Their relative strengths determine the rates of convergence of the estimated factors and loadings.
Under model ((ref)), we can test $\bfsym \gamma_t=0$ almost surely in the entire sampling period, under which the observed $\bfm x_t$ fully explain the true factors. This is the same as testing $$ H_0: \operatorname{cov}(\bfsym \gamma_t)=0. $$ While it is well known that the commonly used Fama-French factors have explanatory power for most of the variations of stock returns, it is questionable whether they fully explain the true (yet unknown) factors. These observed proxies are nevertheless used as the factors empirically, and the remaining components ($\bfsym \gamma_t$ and $\bfm u_t$) have all been mistakenly regarded as the idiosyncratic components. The proposed test provides a diagnostic tool for the specification of common factors in empirical studies, and is different from the “efficiency test" in the financial econometric literature (e.g., GRS, PY, gungor2013testing). While the efficiency test aims to test the asset pricing model through whether the alphas are zero for the specified factors, a rejection could be due to either mispecified factors or the existence of outperforming (underperforming) assets. In contrast, here we directly test whether the factor proxies are correctly specified. We test the specification of Fama French factors for the returns of S&P 500 constituents using rolling windows. We find that the null hypothesis is more often to be rejected using the daily data compared to the monthly data, due to a larger volatility of the unexplained factor components. The estimated overall volatility of factors varies over time and drops significantly during the acceptance period.
In empirical applications, researchers frequently encounter additional observable covariates that help explain the latent factors. In genomic studies, in the study of breast cancer data such as the Cancer Genome Atlas (TCGA) project cancer2012comprehensive, there are additional information of cancer subtype for each sample. These cancer subtypes can be regarded as a partial driver of the factors for gene expression data. In financial time series forecasts, researchers often collect additional variables that characterize financial markets. The Fama-French factors are well-known to be related to the factors that drive financial returns FF.
Most existing works simply treat $\bfm x_t$ as a set of additional regressors in ((ref)). This approach does not take advantage of the difference of observed variables (e.g. aggregated versus disaggregated macroeconomic variables; gene expressions versus clinical information) and the explanatory power of the covariates on the common factors, and hence does not lead to improved rates of convergence even if the signal is strong. The most related work is li2016supervised, who specified $\bfm f_t$ as a linear function of $\bfm x_t$. Also, huang2010combine proposed to use the estimated $\bfm g(\bfm x_t)$ to forecast. Moreover, our expansion is also connected to the literature on asymptotic Bahadur-type representations for robust M-estimators, see, for example, portnoy1985asymptotic, mammen1989asymptotics, among others.
The “asymptotic identification" was described perhaps first by CR. In addition, there has been a large literature on both the static and dynamic factor models, and we refer to Lawley, Forni05, SW02, BN02, bai03, doz, onatski, poet, among many others.
The rest of the paper is organized as follows. Section 2 establishes the new identification of factor models. Section (ref) formally defines our estimators and discusses possible alternatives. Section (ref) presents the rates of convergence. Section (ref) discusses the problem of testing the explanatory power. Section (ref) applies the model to forecasting the excess return of US government bonds. We present the extensive simulation studies in Section (ref) Finally Section (ref) concludes. The supplement also contains all the technical proofs.
Throughout the paper, we use $\lambda_{\min}(\bfm A)$ and $\lambda_{\max}(\bfm A)$ to denote the minimum and maximum eigenvalues of a matrix $\bfm A$. We define $\|\bfm A\|_F=\operatorname{tr}^{1/2}(\bfm A'\bfm A)$, $\|\bfm A\|=\lambda_{\max}^{1/2}(\bfm A'\bfm A)$, $\|\bfm A\|_1 =\max_{j} \sum_{i} |a_{ij}|$ and $\|\bfm A\|_{\max}=\max_{i,j}|a_{ij}|$. For two sequences, we write $a_T\gg b_T$ or $b_T\ll a_T$ if $b_T=o(a_T)$ and $a_T\asymp b_T$ if $a_T=O(b_T)$ and $b_T=O(a_T).$
Suppose that there is a fixed $d$-dimensional observable vector $\bfm x_t$ that is: (i) associated with the latent factors $\bfm f_t$, and (ii) mean-independent of the idiosyncratic term. Taking the conditional mean on both sides of ((ref)), we have
This implies
where
Note that $ E(\bfm y_t|\bfm x_t) $ is identified by the data generating process with observables $\{(\bfm y_t, \bfm x_t)\}_{t\leq T}$, but $\bfsym \Sigma_{f|x}$ is not because $\bfm f_t$ is not observable. Since $N>K$, ((ref)) implies that $\bfsym \Sigma_{y|x}$ is a low-rank matrix, whose rank is at most $K.$ Furthermore, we assume $ \bfsym \Sigma_{f|x}$ is also full rank, so $\bfsym \Sigma_{y|x}$ has exactly $K$ nonzero eigenvalues.
To see how the equality ((ref)) helps achieve the identification of $\ensuremath{\boldsymbol{\Lambda}}$ and $\bfm g(\bfm x_t)$, for the moment, suppose the following normalization holds:
Then right multiplying ((ref)) by $\ensuremath{\boldsymbol{\Lambda}}/N$, by the normalization condition, $$ \frac{1}{N}\bfsym \Sigma_{y|x}\ensuremath{\boldsymbol{\Lambda}} =\ensuremath{\boldsymbol{\Lambda}} \bfsym \Sigma_{f|x}. $$ We see that the ($K$) columns of $\frac{1}{\sqrt{N}}\ensuremath{\boldsymbol{\Lambda}}$ are the eigenvectors of $\bfsym \Sigma_{y|x}$, corresponding to its $K$ nonzero eigenvalues, which also equal to the diagonal entries of $ \bfsym \Sigma_{f|x} $. Furthermore, left multiplying $\ensuremath{\boldsymbol{\Lambda}}'/N$ on both sides of ((ref)), one can see that even if $\bfm f_t$ is not observable, $E(\bfm f_t|\bfm x_t)$ is also identified as: $$ \bfm g(\bfm x_t):=E(\bfm f_t|\bfm x_t)=\frac{1}{N}\ensuremath{\boldsymbol{\Lambda}}' E(\bfm y_t|\bfm x_t). $$
The normalization ((ref)) above is useful to facilitate the above arguments. In this paper, they are not imposed. Then the same argument shows that $\ensuremath{\boldsymbol{\Lambda}}$ and $\bfm g(\bfm x_t)$ can be identified up to a rotation matrix transformation.
Let
Condition (i) is the key condition on the explanatory power of $\bfm x_t$ on factors, where $\chi_N$ represents the “signal strength" of the model. We postpone the discussion of this condition after Theorem (ref). Condition (ii) in Assumption (ref) can be weakened to allow the eigenvalues of $\bfsym \Sigma_{\Lambda, N}$ to slowly decay to zero. While doing so allows some of the factors to be weak, it does not provide any new statistical insights, but would bring unnecessary complications to our results and conditions. Therefore, we maintain the strong version as condition (ii).
Generally, we have the following theorem for identifying $(\ensuremath{\boldsymbol{\Lambda}}, \bfm g(\bfm x_t))$ (up to a rotation transformation).
In the model $$ \bfm f_t=\bfm g(\bfm x_t)+\bfsym \gamma_t,\quad \bfm g(\bfm x_t)=E(\bfm f_t|\bfm x_t), $$ $ \chi_N =\lambda_{\min}(\bfsym \Sigma_{f|x})$ represents the “signal" of the covariate model. We require $\chi_N>0$ so that the rank of $\bfsym \Sigma_{y|x}$ is $K$. Only if this condition holds are we able to identify all the $K$ factor loadings using the eigenvectors corresponding to the nonzero eigenvalues. From the estimation point of view, we are using the PCAs of the estimated $\bfsym \Sigma_{y|x}$, and can only consistently estimate its rank$(\bfsym \Sigma_{y|x})$-number of leading eigenvectors. So this condition is also essential to achieve the consistent estimation of the factor loadings.
Note that requiring $\bfsym \Sigma_{f|x}$ be of full rank might be restrictive in some cases. For instance, consider the linear case: $ E(\bfm f_t|\bfm x_t) =\bfsym \beta\bfm x_t$ for a $K\times d$ coefficient matrix $\bfsym \beta$, also suppose $E\bfm x_t\bfm x_t'$ is of full rank. Then $\bfsym \Sigma_{f|x}=\bfsym \beta E\bfm x_t\bfm x_t' \bfsym \beta'$, and is full-rank only if $d\geq K$. Thus we implicitly require, for linear models, the number of covariates should be at least as many as the number of latent factors. Note that if $E(\bfm f_t|\bfm x_t)$ is highly nonlinear, it is still possible to satisfy the full rank condition even if $d<K$, and we illustrate this in the simulation section. \footnote{Suppose $E(\bfm f_t|\bfm x_t)$ is nonlinear and can be well approximated by a series of orthogonal basis functions $\Phi(\bfm x_t)=(\phi_1(\bfm x_t),...,\phi_J(\bfm x_t))'$, where $E\phi_i(\bfm x_t)\phi_j(\bfm x_t)=1\{i=j\}$, then for some $K\times J$ coefficient $\bfsym \alpha$, we have $E(\bfm f_t|\bfm x_t)\approx \bfsym \alpha' \Phi(\bfm x_t)$ so $E\{E(\bfm f_t|\bfm x_t)E(\bfm f_t|\bfm x_t)'\}\approx \bfsym \alpha\bfsym \alpha'$. For nonlinear functions, it is not stringent to require $\bfsym \alpha\bfsym \alpha'$ be full rank since $K<J$ as $J\to\infty$. }
The above identification strategy motivates us to estimate $\ensuremath{\boldsymbol{\Lambda}}$ and $\bfm g(\bfm x_t)$ respectively by $\widehat{\ensuremath{\boldsymbol{\Lambda}}}$ and $\widehat \bfm g(\bfm x_t)$ as follows. Let $\widehat{\bfsym \Sigma}$ and $\widehat E(\bfm y_t|\bfm x_t)$ be some estimator of $\bfsym \Sigma_{y|x}$ and $E(\bfm y_t|\bfm x_t)$, whose definitions will be clear below. Then the columns of $\frac{1}{\sqrt{N}}\widehat{\ensuremath{\boldsymbol{\Lambda}}}$ are defined as the eigenvectors corresponding to the first $K$ eigenvalues of $\widehat{\bfsym \Sigma}$, and $$\widehat\bfm g(\bfm x_t):=\frac{1}{N}\widehat{\ensuremath{\boldsymbol{\Lambda}}}' \widehat E(\bfm y_t|\bfm x_t). $$ Recall that $\bfm f_t=\bfm g(\bfm x_t)+\bfsym \gamma_t$. We estimate $\bfm f_t$ using least squares: $$ \widehat\bfm f_t:=(\widehat{\ensuremath{\boldsymbol{\Lambda}}}'\widehat{\ensuremath{\boldsymbol{\Lambda}}})^{-1}\widehat{\ensuremath{\boldsymbol{\Lambda}}}'\bfm y_t=\frac{1}{N}\widehat{\ensuremath{\boldsymbol{\Lambda}}}' \bfm y_t. $$ Finally, we estimate $\bfsym \gamma_t$ by: $ \widehat{\bfsym \gamma}_t=\widehat\bfm f_t-\widehat\bfm g(\bfm x_t)=\frac{1}{N}\widehat{\ensuremath{\boldsymbol{\Lambda}}}'(\bfm y_t-\widehat E(\bfm y_t|\bfm x_t)). $ Estimating $\bfm g(\bfm x_t)$ and $\bfsym \gamma_t$ separately allows us to estimate and distinguish the percentage of explained and unexplained components in factors, as well as to quantify the explanatory power of covariates.
Below we introduce the estimators $\widehat{\bfsym \Sigma}$ and $\widehat E(\bfm y_t|\bfm x_t)$ to be used in this paper.
Recall that $\bfsym \Sigma_{y|x}=E\{ E(\bfm y_t|\bfm x_t) E(\bfm y_t|\bfm x_t)' \}$, and let us first construct an estimator for $ E(\bfm y_t|\bfm x_t)$ as follows. While many standard nonparametric regressions would work, here we choose an estimator that is robust to the tail-distributions of $\bfm y_t-E(\bfm y_t|\bfm x_t)$.
Let $\Phi(\bfm x_t)=(\phi_1(\bfm x_t),...,\phi_J(\bfm x_t))'$ be a $J\times 1$ dimensional vector of sieve basis. Suppose $E(\bfm y_t|\bfm x_t)$ can be approximated by a sieve representation: $ E(\bfm y_t|\bfm x_t)\approx \bfm B \Phi(\bfm x_t), $ where $\bfm B=(\bfm b_1,...,\bfm b_N)'$ is an $N\times J$ matrix of sieve coefficients. To adapt to different heaviness of the tails of idiosyncratic components, we use the Huber loss function (huber1964robust) to estimate the sieve coefficients. Define $$ \rho(z)=
$$ For some deterministic sequence $\alpha_T\to\infty$ (adaptive Huber loss), we estimate the sieve coefficients $\bfm B$ by the following convex optimization: $$ \widehat{\bfm b}_i=\arg\min_{b\in\mathbb{R}^J} \frac{1}{T}\sum_{t=1}^T \rho\left( \frac{y_{it}-\Phi(\bfm x_t)'\bfm b}{\alpha_T}\right),\quad \widehat{\bfm b}=(\widehat{\bfm b}_1,...,\widehat{\bfm b}_N)'. $$ We then estimate $\bfsym \Sigma_{y|x}$ by $$ \widehat{\bfsym \Sigma}=\frac{1}{T}\sum_{t=1}^T\widehat E(\bfm y_t|\bfm x_t) \widehat E(\bfm y_t|\bfm x_t)',\quad where \widehat E(\bfm y_t|\bfm x_t)=\widehat{\bfm b}\Phi(\bfm x_t). $$
An alternative method to the robust estimation of $\bfsym \Sigma_{y|x}$ is based on the sieve-least squares, corresponding to the case where $\alpha_T = \infty$. Let $\bfm Y=(\bfm y_1,...,\bfm y_T), $ which is $(N\times T)$, and $$ \bfm P=\Phi'(\Phi\Phi')^{-1}\Phi, (T\times T),\quad \Phi=(\Phi(\bfm x_1),...,\Phi(\bfm x_T)), (J\times T). $$ Then, the sieve least-squares estimator for $\bfsym \Sigma_{y|x}$ is $ \widetilde\bfsym \Sigma= \frac{1}{T} \bfm Y\bfm P\bfm Y'$. While this estimator is attractive due to its closed form, it is not as good as $\widehat\bfsym \Sigma$ when the distribution of $\bfm u_t$ has heavier tails. As expected, our numerical studies in Section (ref) demonstrate that it performs well in light-tailed scenarios, but is less robust to heavy-tailed distributions. Our theories are presented for $\widehat\bfsym \Sigma$, but most of the theoretical findings should carry over to $\widetilde\bfsym \Sigma$.
The selection of the sieve dimension $J$ has been widely studied in the literature, e.g., Li_87, Andrews_91, np_AIC_98, among others. Another tunning parameter is $\alpha_T$, which diverges in order to reduce the biases of estimating the conditional mean when the distribution of $\bfm y_t-E(\bfm y_t|\bfm x_t)$ is asymmetric. Throughout the paper, we shall set
for some constant $C_\alpha>0$, and choose $(J, C_\alpha)$ simultaneously using the multi-fold cross-validation\footnote{One can also allow $\alpha_T$ to depend on $\operatorname{var}(y_{it}|\bfm x_t)$ to allow for different scales across individuals. We describe this choice in the simulation section. In addition, the cross-validation can be based on either in-sample fit for $E(y_{it}|\bfm x_t)$ or out-of-sample forecast, depending on the specific applications. In time series forecasts, one may also consider the time series cross validation TSCV where the training and testing sets are defined through a moving window forecast. }. The specified rate in ((ref)) is due to a theoretical consideration, which leads to the “least biased robust estimation", as we now explain. The Huber-estimator is biased for estimating the mean coefficient in $E(y_{it}|\bfm x_t)$, whose population counterpart is $$ \bfm b_{i,\alpha}:=\arg\min_{\bfm b\in\mathbb{R}^J} E \rho\left( \frac{y_{it}-\Phi(\bfm x_t)'\bfm b}{\alpha_T}\right), $$ As $\alpha_T$ increases, it approaches the limit $\bfm b_i:=\arg\min_{\bfm b\in\mathbb{R}^J} E[y_{it}-\bfm b'\Phi(\bfm x_t)]^2$ with the speed $$ \max_{i\leq N}\|\bfm b_{i,\alpha}-\bfm b_i\|=O(\alpha_T^{-(\zeta_2+1)+\epsilon}) $$ for an arbitrarily small $\epsilon>0$, where $\zeta_2$ is defined in Assumption (ref). Hence the bias decreases as $\alpha_T$ grows. On the other hand, our theory requires the uniform convergence (in $i=1,...,N$) of (for $e_{it}=y_{it}-E(y_{it}|\bfm x_t)$)
where $\dot\rho(\cdot) $ denotes the derivative of $\rho(\cdot)$. It turns out that $\alpha_T$ cannot grow faster than $O(\sqrt{\frac{T}{\log(NJ)}})$ in order to guard for robustness and to have a sharp uniform convergence for ((ref)). Hence the choice ((ref)) leads to the asymptotically least-biased robust estimation.
Plugging $\bfm f_t=\bfm g(\bfm w_t)+\bfsym \gamma_t$ into ((ref)), we obtain
A closely related model is:
for a nonparametric function $\bfm h(\cdot)$, or simply a linear form $\bfm h(\bfm x_t)=\bfsym \beta\bfm x_t$. Models ((ref)) and ((ref)) were studied in the literature ALS, bai, MW11, where parameters are often estimated using least squares. For instance, we can estimate model ((ref)) by
But this approach is not appropriate in the current context when $\bfm x_t$ almost fully explains $\bfm f_t$ for all $t=1,...,T$. In this case, $\bfsym \gamma_t\approx 0$, and least squares ((ref)) would be inconsistent. \footnote{The inconsistency is due to the fact that $a\ensuremath{\boldsymbol{\Lambda}}\bfsym \gamma_t\approx \ensuremath{\boldsymbol{\Lambda}}\bfsym \gamma_t$ for any scalar $a$ in the case $\bfsym \gamma_t\approx0$. Thus $\ensuremath{\boldsymbol{\Lambda}}$ is not identifiable in the least squares problem.} In addition, $\ensuremath{\boldsymbol{\Lambda}}$ in ((ref)) would be very close to zero because the effects of $\bfm f_t$ would be fully explained by $\bfm h(\bfm w_t)$. As a result, the factors in ((ref)) cannot be consistently estimated onatski2012asymptotics either. We conduct numerical comparisons with this method in the simulation section. In all simulated scenarios, the interactive effect approach gives the worst estimation performance.
Another simpler alternative is to combine $(\bfm x_t, \bfm y_t)$, and apply the classical methods on this enlarged dataset. One potential drawback is that the rates of convergence would not be improved, even if $\bfm x_t$ has strong explanatory power on the factors. Another drawback, as mentioned before, is that $\bfm x_t$ and $\bfm y_t$ can provide very different information (e.g. Fama-French factors versus returns of individual stocks).
Let $ e_{it}:=y_{it}-E(y_{it}|\bfm x_t). $ Suppose the conditional distribution of $e_{it}$ given $\bfm x_t=\bfm x$ is absolutely continuous for almost all $\bfm x$, with a conditional density $g_{e, i}(\cdot|\bfm x)$.
Recall $\bfsym \gamma_t=\bfm f_t-E(\bfm f_t|\bfm x_t)$. Let $\gamma_{kt}$ be its $k$ th component.
Recall that $$\bfsym \Sigma_{f|x}:=E\{E(\bfm f_t|\bfm x_t)E(\bfm f_t|\bfm x_t)'\} ,\quad \chi_N:=\lambda_{\min}( \bfsym \Sigma_{f|x} ) .$$
Assumption (ref) allows distributions with relatively heavy tails on $y_{it}-E(y_{it}|\bfm x_t)$. We still require sub-Gaussian tails for the sieve basis functions. Assumption (ref) is regarding the accuracy of sieve approximations for nonparametric functions. Assumption (ref) strengthens Assumption (ref). We respectively regard $ \lambda_{\min}(\bfsym \Sigma_{f|x})$ and $\operatorname{cov}(\bfsym \gamma_t)$ as the “signal" and “noise" when using $\bfm x_t$ to explain common factors. The explanatory power is measured by these two quantities.
Assumption (ref) (i) requires serial independence, and we admit that it can be restrictive in applications. Allowing for serial dependence is technically difficult due to the non-smooth Huber's loss. To obtan the Bahadur representation of the estimated eigenvectors, we rely on the symmetrization and contraction theorems (e.g., VW), which requires the data be independently distributed. Nevertheless, the idea of using covariates would still be applicable for serial dependent data. For instance, it is not difficult to allow for weak serial correlations when the data are not heavy-tailed, by using the sieve least squares estimator $\widetilde\bfsym \Sigma$ (introduced in Section (ref)) in place of the Huber's estimator $\widehat{\bfsym \Sigma}$. We conduct numerical studies when the data are serially correlated in the simulations, and find that the proposed methods continue to perform well in the presence of serial correlations.
We present the rates of convergence in the following theorems, and discuss the statistical insights in the next subsection. Recall $\widehat{\ensuremath{\boldsymbol{\Lambda}}}=(\widehat{\ensuremath{\boldsymbol{\lambda}}}_i: i\leq N)$.
The optimal rate for $J$ in ((ref)) is $J\asymp T^{1/(2\eta)}$, which results in
Here $\eta$ represents the smoothness of $E(\bfm f_t|\bfm x_t=\cdot)$, as defined in Assumption (ref).
Define $$ J^{*}= \min\left\{ (TN)^{1/(2\eta)}, ( \frac{T}{\log N})^{1/(1+\eta)} \right\}. $$
These two convergences imply the rate of convergence of the estimated factors due to $\widehat\bfm f_t=\widehat \bfm g(\bfm x_t)+\widehat{\bfsym \gamma}_t$.
We see that the rates depend on $\operatorname{cov}(\bfsym \gamma_t)$ and $\chi_N$. Because $E\bfm f_t\bfm f_t'=\bfsym \Sigma_{f|x}+\operatorname{cov}(\bfsym \gamma_t)$, they are related through
for some $c, C_1 >0$, assuming that there is $c>0$ so that $\|E\bfm f_t\bfm f_t'\|>c$. For comparison, we state the rates of convergence of the benchmark PCA estimators: (e.g., SW02, bai03) there is a rotation matrix $\tilde\bfm H$, so that the PCA estimators $( \widetilde\ensuremath{\boldsymbol{\lambda}}_i, \widetilde\bfm f_t)$ satisfy:
The first interesting phenomena we observe is that both the estimated loadings and $\bfm g(\bfm x_t)$ are consistent even if $N$ is finite, due to the “exact identification". In contrast, the PCA estimators requires a growing $N$. For more detailed comparisons, we consider three regimes based on the explanatory power of the factors using $\bfm x_t$. To simplify our discussions, we consider the rate-optimal choices of $J$, and ignore the sieve approximation errors, so $\eta$ is treated sufficiently large.
Regime I: strong explanatory power: $\|\operatorname{cov}(\bfsym \gamma_t)\|\to 0.$ Because of ((ref)), $\chi_N$ is bounded away from zero. In this case, ((ref))-((ref)) approximately imply (for sufficiently large $\eta$):
Compared to the rates of the usual PCA estimators in ((ref)), either the new estimated loadings (when $N=o(T)$) or the new estimated factors (when $T=o(N)$) have a faster rate of convergence. Moreover, if $\|\operatorname{cov}(\bfsym \gamma_t)\|= o((TN)^{-1}+T^{-2}\log^2N)$, then $\widehat\bfm g(\bfm x_t)$ directly estimates the latent factor at a very fast rate of convergence: $$ \frac{1}{T}\sum_{t=1}^T\|\widehat\bfm g(\bfm x_t)-\bfm H^{-1}\bfm f_t\|^2=O_P\left( \frac{1}{TN}+ (\frac{\log N}{T})^{2 } \right). $$ The improved rates are reasonable due to the strong explanatory powers from the covariates.
Regime II: mild explanatory power: $\|\operatorname{cov}(\bfsym \gamma_t)\| $ is bounded away from zero; $\chi_N$ is either bounded away from zero or decays slower than $ \frac{N}{T}$ in the case $N=o(T)$. In this regime, $\bfm x_t$ partially explains the factors, yet the unexplainable components are not negligible. ((ref))-((ref)) approximately become:
We see that the rate for the estimated loadings is still faster than the PCA when $N$ is relatively small compared to $T$, while the rates for the estimated factors are the same. This is because $$ \underbrace{\frac{1}{T}\chi_N^{-1}}_{\text{new rate for loadings}}\ll \underbrace{\frac{1}{N}}_{\text{ PCA rate for loadings}} $$ $$ \underbrace{\frac{1}{T}\chi_N^{-1}+ \frac{1}{N}}_{\text{new rate for factors}}\asymp \underbrace{\frac{1}{N}}_{\text{ PCA rate for factors}}. $$ On one hand, due to the explanatory power from the covariates, the loadings can be estimated well without having to consistently estimate the factors. On the other hand, as the covariates only partially explain the factors, we cannot improve rates of convergence in estimating the unexplainable components in the latent factors. However, since $\bfsym \gamma_t$ has smaller variability than $\bfm f_t$, it can still be better estimated in terms of a smaller constant factor.
Regime III: weak explanatory power: $\chi_N\to 0$ and decays faster than $ \frac{N}{T}$ when $N\ll T$. In this case, we have
While the new estimators are still consistent, they perform worse than PCA. This finding is still reasonable because the signal is so weak that the conditional expectation $E(\bfm y_t|\bfm x_t)$ loses useful information of the factors/loadings. Consequently, estimation efficiency is lost when running PCA on the estimated covariance $ E\{ E(\bfm y_t|\bfm x_t) E(\bfm y_t|\bfm x_t)' \}$.
In summary, improved rates of convergence can be achieved so long as the covariates can (partially) explain the latent factors, this corresponds to either the mild or the strong explanatory power case. The degree of improvements depend on the strength of the signals. In particular, the consistent estimation for factor loadings can also be achieved even under finite $N$. On the other hand, when the explanatory power is too weak, the rates of convergence would be slower than those of the benchmark estimator.
We aim to test: (recall that $\bfsym \gamma_t=\bfm f_t-E(\bfm f_t|\bfm x_t)$)
Under $H_0$, $\bfm f_t=E(\bfm f_t|\bfm x_t)$ over the entire sampling period $t=1,..., T$, implying that observed covariates $\bfm x_t$ fully explain the true factors $\bfm f_t$. In empirical applications with “observed factors", what have been often used are in fact $\bfm x_t$. Hence our proposed test can be applied to empirically validate the explanatory power of these “observed factors".
The Fama-French three-factor model FF is one of the most celebrated ones in empirical asset pricing. They modeled the excess return $r_{it}$ on security or portfolio $i$ for period $t$ as $$ r_{it} =\alpha_i +b_i r_{Mt} + s_i \mbox{SMB}_t +h_i \mbox{HML}_t + u_{it}, $$ where $r_{Mt}, \mbox{SMB}_t$ and $\mbox{HML}_t$ respectively represent the the excess returns of the market, the difference of returns between stocks with small and big market capitalizations (“small minus big"), and the difference of returns between stocks with high book to equity ratios and those with low book to equity ratios (“high minus low"). Ever since its proposal, there is much evidence that the three-factor model can leave the cross-section of expected stock returns unexplained. Different factor definitions have been explored, e.g., carhart1997persistence and novy2013other. fama2015five added profitability and investment factors to the three-factor model. They conducted GRS tests GRS on the five-factor models and its different variations. Their tests “reject all models as a complete description of expected returns".
On the other hand, the Fama-French factors, though imperfect, are good proxies for the true unknown factors. Consequently, they form a natural choice for $\bfm x_t$. These observables are actually diversified portfolios, which have explanatory power on the latent factors $\bfm f_t$, as supported by financial economic theories as well as empirical studies. The test proposed in this validates the specification of these common covariates as “factors".
Our test is based on a Wald-type weighted quadratic statistic $$ S(\bfm W):=\frac{N}{T}\sum_{t=1}^T\widehat{\bfsym \gamma}_t'\bfm W\widehat{\bfsym \gamma}_t=\frac{1}{TN}\sum_{t=1}^T(\bfm y_t-\widehat E(\bfm y_t|\bfm x_t))'\widehat{\ensuremath{\boldsymbol{\Lambda}}}\bfm W\widehat{\ensuremath{\boldsymbol{\Lambda}}}'(\bfm y_t-\widehat E(\bfm y_t|\bfm x_t)). $$ The weight matrix normalizes the test statistic, taken as $\bfm W=$ AVar$(\sqrt{N}\widehat{\bfsym \gamma}_t)^{-1}$, where AVar$(\widehat{\bfsym \gamma}_t)$ represents the asymptotic covariance matrix of $\widehat{\bfsym \gamma}_t$ under the null, and is given by $$ \text{AVar}(\sqrt{N}\widehat{\bfsym \gamma}_t)=\frac{1}{N}\bfm H' \ensuremath{\boldsymbol{\Lambda}}'\bfsym \Sigma_u\ensuremath{\boldsymbol{\Lambda}}\bfm H. $$ As $\bfsym \Sigma_u$ is a high-dimensional covariance matrix, to simplify the technical arguments, in this section we assume $\{u_{it}\}$ to be cross-sectionally uncorrelated, and estimate $\bfsym \Sigma_u$ by: $$ \widehat{\bfsym \Sigma}_u=\operatorname{diag}\{\frac{1}{T}\sum_{t=1}^T\widehat u_{it}^2, i=1,..., N\}, \quad \widehat u_{it}=y_{it}-\widehat{\ensuremath{\boldsymbol{\lambda}}}_i'\widehat\bfm f_t. $$ The feasible test statistic is defined as $$ S := S(\widehat\bfm W),\quad \widehat\bfm W:=(\frac{1}{N}\widehat{\ensuremath{\boldsymbol{\Lambda}}}'\widehat{\bfsym \Sigma}_u\widehat{\ensuremath{\boldsymbol{\Lambda}}})^{-1}. $$ We reject the null hypothesis for large values of $S$. It is straightforward to allow $\bfsym \Sigma_u$ to be a non-diagonal but a sparse covariance, and proceed as in Bickel08a. We expect the asymptotic analysis to be quite involved, and do not pursue it in this paper.
We show that the test statistic has the following asymptotic expansion: $$ S= \bar S+ o_P(\frac{1}{\sqrt{T}}), $$ where $$ \bar S=\frac{1}{T}\sum_{t=1}^T \bfm u_t'\ensuremath{\boldsymbol{\Lambda}} (\ensuremath{\boldsymbol{\Lambda}}'\bfsym \Sigma_u\ensuremath{\boldsymbol{\Lambda}})^{-1}\ensuremath{\boldsymbol{\Lambda}}'\bfm u_t. $$ Thus the limiting distribution is determined by that of $\bar S$. Note that a cross-sectional central limit theorem implies, as $N\to\infty$, $$ (\frac{1}{N}\ensuremath{\boldsymbol{\Lambda}}'\bfsym \Sigma_u\ensuremath{\boldsymbol{\Lambda}})^{-1/2}\frac{1}{\sqrt{N}}\bfm u_t'\ensuremath{\boldsymbol{\Lambda}}\to^d\mathcal{N}(0,\bfm I_K). $$ Hence each component of $\bar S$ can be roughly understood as $\chi^2$-distributed with degrees of freedom $K$ being the number of common factors, whose variance is $2K$. This motivates the following assumption.
We now state the null distribution in the following theorem.
\numberwithin{table}{section} \numberwithin{figure}{section}
We test the explanatory power of the observable proxies for the true factors using S&P 500 returns. For each given group of observable proxies, we set the number of common factors $K$ equals the number of observable proxies. We calculate the excess returns for the stocks in S&P 500 index that are collected from CRSP. We consider three groups of proxy factors ($\bfm x_t$) with increasing information: (1) Fama-French 3 factors (FF3); (2) Fama-French 5 factors (FF5); and (3) Fama-French 5 factors plus 9 sector SPDR ETF's (FF5+ETF9). Here the sector SPDR ETF's, which are intended to track the 9 largest S&P sectors. The detailed descriptions of sector SPDR ETF's are listed in Table (ref).
We consider tests using both daily and monthly data. For the daily data, we collect 393 stocks that have complete daily closing prices from January 2005 to December 2013, with a time span of 2265 trading days. We apply moving window tests with the window size ($T$) equals one month, three months or six months. The testing window moves one trading day forward per test. Within each testing window, we calculate the standardized test statistic $S$ for three groups of proxy factors.
As for the monthly excess returns, we use stocks that have complete record from January 1980 to December 2012, which contains 202 stocks with a time span of 396 months. Here we only consider the first two groups of proxy factors as sector SPDR ETF's are introduced since 1998. The window size equals sixty months and moves one month forward per test. Within each testing window, besides standardized test statistic and p-value, we also estimate the volatility of $\bfsym \gamma_t$, the part of factors that can not be explained by $\bfm x_t$ as: $$ \widehat{\mbox{Vol}}(\bfsym \gamma_t)=\frac{1}{21T}\sum\limits_{t=1}^T\widehat\bfsym \gamma_t'\widehat\bfsym \gamma_t, $$ where there are 21 trading days per month. The sieve basis is chosen as the additive Fourier basis with $J=5$. We set the tuning parameter $\alpha_T=C\sqrt{\frac{T}{\log(NJ)}}$ with constant $C$ selected by the 5-fold cross validation.
For the daily data, the plots of $S$ under various scenarios are reported in Figure (ref). Under all scenarios, the null hypothesis ($H_0: \operatorname{cov}(\bfsym \gamma_t)=0$) is rejected as $S $ is always larger than the critical value 1.96. This suggests a strong evidence that the proxy factors can not fully explain the estimated common factors. Under all window sizes, a larger group of proxy factors tends to yield smaller statistics, demonstrating stronger explanatory power for estimated common factors. Also, we find the test statistics increase while the window size increases.
The results for the monthly data are reported in Figure (ref). For both Fama-French 3 factors and 5 factors, the null hypothesis is rejected most of the time except in early 1980s and 1990s. When the null hypothesis is accepted, Fama-French 5 factors tend to yield larger p-values. The estimated volatility of unexplained part are close to zero over these two periods. For the rest of the time, the standardized test statistics are much larger than the critical value 1.96 and hence the p-values are close to zero. Also the estimated volatilities are not close to zero. This indicates the proxy factors can not fully explain the estimated common factors during these testing periods.
We apply our method to forecast the excess return of U.S. government bonds. The bond excess return is the one-year bond return in excess of the risk-free rate. To be more specific, we buy an $n$ year bond, sell it as an $n-1$ year bond in the next year and excess the one-year bond yield as the risk-free rate. Let $p_t^{(n)}$ be the log price of an $n$-year discount bond at time $t$. Denote $\zeta_t^{(n)} \equiv -\frac{1}{n}p_t^{(n)}$ as the log yield with $n$ year maturity, and $r_{t+1}^{(n)} \equiv p_{t+1}^{(n-1)}-p_t^{(n)}$ as the log holding period return. The goal of one-step-ahead forecast is to forecast $z_{T+1}^{(n)}$, the excess return with maturity of $n$ years in period $T+1$, where $$ z_{t+1}^{(n)}=r_{t+1}^{(n)}-\zeta_t^{(1)},\quad t=1, \ \cdots \ , T.$$
For a long time, the literature has found a significant predictive power of the excess returns of U.S. government bonds. For instance, ludvigson2009macro, ludvigson2010factor predicted the bond excess returns with observable variables based on a factor model using 131 (disaggregated) macroeconomics variables. They achieved the out-of-sample $R^2 \approx 21\%$ when forecasting one year excess bond return with maturity of two years. Using the proposed method, this section develops a new way of incorporating the explanatory power of the observed characteristics, and investigates the robustness of the conclusions in existing literature.
We analyze monthly data spanned from January 1964 to December 2003, which is available from the Center for Research in Securities Prices (CRSP). The factors are estimated from a macroeconomic dataset consisting of 131 disaggregated macroeconomic time series ludvigson2010factor. The covariates $\bfm x_t$ are 8 aggregated macro-economic time series, listed in Table (ref).
We first examine the excess kurtosis for the time series to assess the tail distributions. The left panel of Figure (ref) shows 43 among the 131 series have excess kurtosis greater than 6. This indicates the tails of their distributions are fatter than the $t$-distribution with degrees of freedom 5. On the other hand, the right panel of Figure (ref) reports the histograms of excess kurtosis of the “fitted data" $\widehat{E}(\bfm y_t|\bfm x_t)$ (the robust estimator of $E(\bfm y_t|\bfm x_t)$ using Huber loss), which demonstrates that most series in the fitted data are no longer severely heavy-tailed.
The tuning parameter in the Huber loss is of order $\alpha_{T}=C_\alpha \sqrt{\frac{T}{\log(NT)}}$. In this study, the constant $C_{\alpha}$ and the degree of sieve approximation $J$ are selected by the out-of-sample 5-fold cross validation as described in Section (ref).
We denote our proposed method by SPCA (smoothed PCA), and compare it with SPCA-LS (which uses $\widetilde\bfsym \Sigma$, the least-squares based smoothed PCA, described in Section (ref)) and the benchmark PCA. We conduct one-month-ahead out-of-sample forecast of the bond excess returns. The forecast uses the information in the past 240 months, starting from January 1984 and rolling forward to December 2003. We compare three approaches to estimating the factors: SPCA, SPCA-LS, and the usual PCA. Also we consider two forecast models as follows:
where $\alpha$ is the intercept and $h$ is a nonparametric function. The covariate $\bfm W_t$ is either $\bfm f_t$ or an augmented vector $(\bfm f_t',\bfm x_t')'$. Here, the latent factors $\bfm f_t$ are used by the three methods mentioned above in order to compare their effectiveness. The multi-index model allows more general nonlinear forecasts and are estimated by using the sliced inverse regression li1991sliced. The number of indices $L$ is estimated by the ratio-based method suggested in LamYao and is usually 2 or 3. We approximate $h$ using a weighted additive model $h({\ensuremath{\boldsymbol{\psi}}}_1' \bfm W_t,\ \cdots , \ {\ensuremath{\boldsymbol{\psi}}}_L' \bfm W_t)=\sum_{l=1}^L g_l({\ensuremath{\boldsymbol{\psi}}}_l' \bfm W_t)$. Each individual nonparametric function $g_l(\cdot)$ is smoothed by the local linear approximation.
The performance of each method is assessed by the out-of-sample $R^2$. Let $\widehat z_{T+t+1|T+t}$ be the forecast of $z_{T+t+1}$ using the data of the previous $T$ months: $1+t,...,T+t$ for $T=240$ and $t=0,..., 239$. The forecast performance is assessed by the out-of-sample $R^2$, defined as $$ R^2=1- \frac{\sum\limits_{t=0}^{239}(z_{T+t+1}-\widehat{z}_{T+t+1|T+t})^2}{\sum\limits_{t=0}^{239}(z_{T+t+1}-\bar{z}_t)^2}, $$ where $\bar{z}_t$ is the sample mean of $z_t$ over the sample period $[1+t, T+t]$. The $R^2$ of various methods are reported in Table (ref). We notice that factors estimated by SPCA and SPCA-LS can explain more variations in bond excess returns with all maturities than the ones estimated by PCA. SPCA yields a 44.6% out-of-sample $R^2$ for forecasting the bond excess returns with two year maturity, which is much higher than the best out-of-sample predictor found in ludvigson2009macro. It is also observed that the forecast based on either SPCA or SPCA-LS cannot be improved by adding any covariate in $\bfm x_t$. We argue that, in this application, the information of $\bfm x_t$ should be mainly used as the explanatory power for the factors.
We summarize the observed results in the following aspects:
We use simulated examples to demonstrate the finite sample performance of the proposed method, which is denoted by SPCA (smoothed PCA), and compare it with SPCA-LS (which uses $\widetilde\bfsym \Sigma$, the least-squares based smoothed PCA, described in Section (ref)) and the benchmark PCA. We set $N=40, \ T=100$ and $K=5$. The supplementary material contains additional simulation results under other $N$ and $T$ combinations, as well as the case of serially dependent data. The findings are similar.
Consider the following data generating process,
where $\ensuremath{\boldsymbol{\Lambda}}$ is drawn from i.i.d. standard Normal distribution and $\bfm u_t$ is drawn from either the i.i.d standard Normal distribution or i.i.d. re-scaled Log-Normal distribution $c_1 \{ \exp(1+1.2 \zeta ) - c_2\}$, where $\zeta \sim \mathcal{N}(0, 1)$ and $c_1, c_2 >0$ are chosen such that $u_{it}$ has mean zero and variance $1$.
Here $\tilde{\sigma}(g)$ and $\tilde{\sigma}(\gamma)$ respectively represent the signal and noise levels. Set $\tilde{\sigma}(g)^2+\tilde{\sigma}(\gamma)^2=1$ and $\tilde{\sigma}(g)^2/\tilde{\sigma}(\gamma)^2=\omega$, where $\omega$ controls the ratio between the explained and unexplained parts in the latent factors. To address different signal-noise regimes, we set $\omega=10, 1$ and $0.1$ to represent strong, mild and weak explanatory powers respectively. The baseline $\bfsym \gamma_t^0$ is drawn from i.i.d. standard Normal distribution and the baseline function $\bfm g^0(\cdot)$ is set to be one of the following two models:
For each $k\leq K$, we normalize $g_k^0(\bfm x_t)$ and $\gamma^0_{t,k}$ such that they have means zero, and standard deviations one.
Throughout this section, the number of factors is estimated by the eigen-ratio method LamYao,AH. In the following simulated examples, the eigen-ratio method can correctly select $K=5$ in most replications. The sieve basis is chosen as the additive polynomial basis. To account the scale of the noise variance, we also consider the tuning parameter in the Huber loss to admit $\alpha_{T, i}=C_\alpha\tilde{\sigma}_i\sqrt{\frac{T}{\log(NT)}}$, where $\tilde{\sigma}_i=\sqrt{ \frac{1}{T}\sum_t(y_{it}-\tilde E(y_{it}|\bfm x_t))^2}$ and $\tilde E(y_{it}|\bfm x_t)$ is smoothed by sieve least squares using additive polynomial basis of order 5. In Subsection (ref), the tuning parameters $C_{\alpha}$ and $J$ are selected by the in-sample 5-fold cross validation, while in subsection (ref), they are chosen using the out-of-sample 5-fold cross validation.
First, we compare the in-sample model fitting among SPCA, SPCA-LS and PCA under different scenarios. For each scenario, we conduct 200 replications. As the factors and loading may be estimated up to a rotation matrix, the canonical correlations between the parameter and its estimator can be used to measure the estimation accuracy bai03. For Model (I) and (II) we report the sample mean of the median of 5 canonical correlations between the true loading and factors and the estimated ones.
The results are presented in Table (ref). SPCA-LS and SPCA are comparable for light-tail distributions, and are both slightly better than PCA. This implies that we pay little price for the robustness and that the proposed estimators are potentially better than PCA when $N$ is relatively small, due to the merit of the “finite-$N$” asymptotics of the proposed estimators. However, when the error distributions have heavy tails, SPCA yields much better estimation than other methods as expected. SPCA-LS out-performs PCA when $\bfm x_t$ has strong or mild explanatory powers of $\bfm f_t$ which is in line with the discussion in Section (ref). When $\omega = 0.1$, the observed $\bfm x_t$ is not as informative and hence the performance of SPCA and SPCA-LS are close to regular PCA.
We now consider using latent factors in a linear forecast model $ z_{t+1} =\bfsym \beta' \bfm f_t + \epsilon_{t+1}, $ where $\epsilon_{t}$ is drawn from i.i.d. standard normal distribution. For each simulation, the unknown coefficients in $\bfsym \beta$ are independently drawn from uniform distribution $[0.5, \ 1.5]$ to cover a variety of model settings.
We conduct one-step ahead rolling window forecast using the linear model by estimating $\bfsym \beta$ and $\bfm f_t$. The factors are estimated from ((ref)) by SPCA, SPCA-LS or PCA. In each replication, we generate $T+50$ observations in total. For $s=1, \ \cdots \ ,50$, we use the $T$ observations $(z_s,..., z_{T+s-1})$ to forecast $z_{T+s}$. We use PCA as the benchmark and define the relative mean squared error (RMSE) as:
where $\widehat{z}_{T+s|T+s-1} $ is the forecast of $z_{T+s}$ based on either SPCA or SPCA-LS while $\tilde{z}_{T+s|T+s-1}^{PCA}$ is the forecast based on PCA.
For each scenario, we simulate 200 replications and calculate the averaged RMSE as a measurement of the one-step-ahead out-of-sample forecast.
The results are presented in Table (ref). Again, when the tails of error distributions are light, SPCA and SPCA-LS perform comparably. But SPCA outperforms SPCA-LS when the errors have heavy tails. On the other hand, both SPCA and SPCA-LS outperform PCA when $\bfm x_t$ exhibits strong or mild explanatory powers of $\bfm f_t$, but are slightly worse when $\omega$ is small. In general, the SPCA method performs the best under heavy-tailed cases.
Here we consider three pairs of sample sizes: $N=40,\ T=150$; $N=60,\ T=100$ and $N=60, \ T=150$. We compare the proposed SPCA method with SPCA-LS (Section (ref)), regular PCA and pure least squares (LS), which models the covariates and estimates the parameters by simply using $$ \min_{\ensuremath{\boldsymbol{\Lambda}},\{\bfm f_t\},\bfsym \beta}\frac{1}{T}\sum_{t=1}^T\|\bfm y_t-\ensuremath{\boldsymbol{\Lambda}} \bfm f_t-\bfm x_t'\bfsym \beta\|^2. $$
In Tables (ref)--(ref), we report sample mean of the median of 5 canonical correlations between the true loading and factors and the estimated ones. Under various sample size combinations, the findings are similar as discussed in Section (ref): (1) both SPCA and SPCA-LS outperform PCA under light-tail distributions when $\bfm x_t$ has strong or mild explanatory powers of $\bfm f_t$; (2) when the error distributions have heavy tails, SPCA outperforms other methods as expected; (3)when $\bfm x_t$ has weak explanatory power, the performance of SPCA and SPCA-LS are close to regular PCA; (4) under all simulated scenarios, the LS approach gives the worst estimation performance.
In this subsection, we compare the in-sample model fitting among SPCA, SPCA-LS and PCA under serial dependences. The simulation settings are similar as in Section 5.1 except both $\bfm x_t$ and $\bfsym \gamma_t$ are generated from a stationary VAR(1) model as follows
with $\bfm x_0={\bf 0}$ and $\bfsym \gamma_0= {\bf 0}$. The $(i, j)$th entry of ${\bf \Pi}$ is set to be 0.5 when $i=j$ and $0.4^{|i-j|}$ when $i \neq j$. In addition, $\bfsym \varepsilon_t$ and ${\bf \eta}_t$ are drawn form i.i.d. $N({\bf 0}, \bfm I)$.
The performance under 200 replications are presented in Table (ref) below. Our numerical findings for the independent data continue to hold for serially dependent data: both SPCA and SPCA-LS outperform PCA when $\bfm x_t$ and $\bfm f_t$ are serially correlated. SPCA gives the best performance when the error distributions are heavy-tailed.
We study factor models when the factors depend on observed explanatory characteristics. The proposed method incorporates the explanatory power of these observed covariates, and is robust to possibly heavy-tailed distributions. We focus on the case $\dim(\bfm x_t)$ is finite, and on the rates of convergence for the estimated factors and loadings. Under various signal-noise ratios, substantial improved rates of convergence can be gained.
Related to the above, the idea could be easily extended to the case that $\dim(\bfm x_t)$ is slowly growing (with respect to $(N, T)$). On the other hand, allowing $\dim(\bfm x_t)$ to be fast-growing would require some dimension-reduction treatment combined with covariate selections. In addition, selecting the covariates would be also useful as the quality of the signal is crucial. We shall leave these open questions for future studies.