Forward Selection Fama-MacBeth Regression with Higher-Order Asset Pricing Factors

Higher-Order Asset Pricing Factors via Forward Selection Fama-MacBeth Regression

\doublespacing

\noindentKeywords: high-dimensional factor models, risk premium, stochastic discount factor, forward selection Fama-MacBeth regression procedure, factor zoo

Introduction

The asset pricing literature has proposed hundreds of candidate factors to explain the cross-section of asset returns. cochrane2011presidential, in his presidential address, refers to this proliferation of candidate factors as the “factor zoo.” This proliferation is itself symptomatic of the fact that sparse linear factor models, such as those introduced by fama2015five and Carhart1997, struggle to fully price the cross-section and hence may not fully span the stochastic discount factor (SDF). In this paper, we show that the interactions of the higher-order terms of the common sparse linear factors are significantly priced in the cross-section of equity returns and capture a significant fraction of the factor zoo. We do so by determining the important interactions via a new factor-selection procedure, which we call the forward selection Fama-MacBeth regression procedure (FS-FMB).

The prevalent focus in the literature on sparse linear factor models is primarily driven by considerations of parsimony and is motivated by the first-order Taylor expansion of the SDF. This observation naturally raises the question of whether higher-order terms and interactions of common sparse linear factors, motivated by a higher-order Taylor expansion, also play an important role in pricing the cross-section of expected returns. In this paper, we show that a small and interpretable subset does, by selecting and ranking the higher-order terms.

A main challenge in studying such relevant terms is dimensionality. Consider, for example, the standard 6-factor model that includes the Fama-French 5 factors (fama2015five, fama2015five) and the momentum factor (Carhart1997, Carhart1997). The higher-order interactions up to degree 3 (4) already have 63 (114) factors. Therefore, we develop a forward selection procedure to choose factors that explain the stochastic discount factor that is suitable for the case when the number of potentially important factors is large. The method adds factors one by one into the set of selected factors, each time maximizing the $R^{2}$ of the cross-sectional regression (second step) of the Fama-MacBeth procedure, thus naturally ranking the importance of the selected factors. After the set of factors is selected, we estimate the SDF loadings by running the cross-sectional regression of sample average returns on sample covariances between returns and selected factors. We derive the rate of convergence of the proposed estimator and also propose a method to debias it, leading to simple inference procedures.

We show that the performance of the common factors is substantially improved with their higher-order interactions. We start with 6 of the most common factors in the asset pricing literature. These are the Fama-French 5 factors of fama2015five and the momentum factor of jegadeesh1993returns, and we refer to this model as the FF5M. We consider higher-orders as well as interactions of these 6 factors up to degree 3, leading to 57 candidate higher-order factors. For the test assets, we use a large cross-section of equities including 484 characteristic-managed portfolios from kozak2020shrinking. We apply our FS-FMB to select the higher-order factors that are most important in pricing the cross-section of these 484 asset returns. The procedure selects 7 higher-order factors, including 2 second degree higher-orders and 5 interactions.\footnote{The FS-FMB procedure requires the choice of a stopping threshold for the second-stage $R^2$ improvement. In our baseline, we set it to be 1%. In the Appendix, we show that we can use a cross-validation method to choose the optimal stopping threshold, which generates qualitatively and quantitatively similar results.} The inclusion of the 5 interaction terms shows the importance of co-dependence of common asset pricing factors in pricing the cross-section of equity returns. We show that the inclusion of the higher-order factors leads to substantial improvement in pricing the cross-section of asset returns -- the adjusted cross-sectional $R^{2}$ increases from 31.2% for the baseline FF5M model to 59% for the selected higher-order factor model. Moreover, the intercept alpha of the higher-order factor model is not statistically significant at the 5% level.

We then examine the SDF loadings of the higher-order factors selected by the FS-FMB using both the non-debiased estimates based on the standard Fama-MacBeth (FMB) two-pass regression and the debiased estimates. Estimates from the standard FMB two-pass procedure are likely biased because of the omitted variable bias (see, e.g., chernozhukov2018double and giglio2021asset). The directions of the point estimates are the same and the magnitudes are similar across the two procedures. The estimates based on the non-debiased procedure are all significant at the 5% level, while the estimates based on the debiased procedure are all significant at the 10% level with 5 of the 7 SDF loading estimates being significant at the 5% level. Overall, the results show that the higher-order factors are important for empirically capturing the SDF and significantly price assets in the cross-section.

In our baseline analysis, we demonstrate that higher-order factors are important components of the SDF and price cross-sectional asset returns in-sample. We further provide evidence of the out-of-sample performance of the selected higher-order factor model, where we conduct two types of out-of-sample tests. Firstly, we present an out-of-sample test in the asset space. We select higher-order factors and estimate SDF loadings from the training assets and then evaluate the out-of-sample performance using the remaining assets. Our results show that the selected higher-order factors and the performance of the higher-order factor model are stable in this out-of-sample exercise.

Secondly, we conduct an out-of-sample test in the time series, where we split the sample into two equal subsamples based on the time series. We use the first half of the sample as the training sample and the second half of the sample as the test sample. Then, we estimate the model in the training sample and evaluate the out-of-sample performance in the test sample. Compared with the benchmark factor models including the CAPM, the FF3 factor model, the FF5 factor model, and the FF5M factor model, the out-of-sample $R^{2}$ of the higher-order factor model is significantly larger.

Next, we examine whether the selected higher-order factor model can reduce the pricing power of the factors in the factor zoo. We employ two empirical procedures. First, we perform standard FMB regressions to test the significance of the SDF loading of each factor in the factor zoo, controlling for the higher-order factor model. Second, we construct factor mimicking portfolios for the seven higher-order factors and conduct time-series asset pricing tests on the zoo factors.

Under the first approach, we estimate each zoo factor's SDF loading while controlling for our higher-order factor model. We find that only 7 out of the 148 zoo factors from jensen2023there remain significant at the 5% level. The second approach is based on factor-mimicking portfolios. The benefit is that we are able to convert the non-tradable higher-order factors to tradable factors but the projection inevitably loses information from the original factors. The results using the second approach confirm that the higher-order factor model significantly reduces the pricing power of the zoo factors. Overall, these results demonstrate that the higher-order factors selected by the forward selection Fama-MacBeth procedure substantially reduce the zoo factors, although the linear benchmark factors cannot.

Because our FS-FMB procedure directly selects and ranks the important higher-order factors while preserving their interpretability, one of its key advantages is that it allows us to connect these factors to underlying economic mechanisms. To this end, we first investigate the economic content of the selected higher-order factors by examining their exposures to a broad cross-section of macroeconomic variables. While unconditional regressions reveal some baseline sensitivity, especially to financial uncertainty, we find that the economic content of these factors is concentrated in extreme economic states. Using quantile regressions, we show that both the explanatory power and the sensitivity to specific macroeconomic variables are markedly amplified in the tails of the distribution. For example, the pseudo-$R^2$ for the market factor squared reaches 70% at the 90th percentile, and variables related to intermediary capital constraints, which are often insignificant in unconditional regressions, become strongly significant in the tails. These results suggest that the higher-order factors capture non-linear exposure to macroeconomic tail risks, particularly those associated with financial intermediary distress. This provides a potential economic rationale, related to intermediary asset pricing (e.g., he2013intermediary,brunnermeier2014macroeconomic), for why higher-order terms of common factors are priced and why they account for a substantial fraction of the broader factor zoo.

Second, we study how these higher-order factors vary with broader macroeconomic and business-cycle conditions by examining their correlations with two key indicators, U.S. equity market returns and the CBOE VIX index, both in the full sample and during NBER recession periods. For the first two higher-order factors (i.e., SMB2 and SMB2{*}Mom), the magnitude of the correlations with U.S. equity market returns is substantially larger during recessions than in the full sample. For the last four higher-order factors (i.e., Mkt-RF2, Mkt-RF2{*}RMW, Mkt-RF{*}SMB, and HML2{*}Mkt-RF), the correlations with the VIX index are large in magnitude, especially during recessions. Overall, these results show that the selected higher-order factors are closely linked to macroeconomic risk over the business cycle, with particularly strong exposure during downturns.

Additionally, we provide robustness results for our empirical analyses. First, we show that the performance of the higher-order factor model is extremely unlikely to be driven by randomness by simulating white-noise factors. Second, we use the selected higher-order factor model to account for an alternative factor zoo of ChenZimmermann2021 and reach similar conclusions. Third, we experiment with alternative sets of candidate higher-order factors: (1) include only higher-orders, (2) only interactions, (3) up to degree 2, and (4) up to degree 4. We find that both higher-orders and interactions are important components of the SDF, especially interactions. Moreover, we show that the selected higher-order factor model up to degree 2 significantly underperforms, while the selected higher-order factor model up to degree 4 has similar performance to our baseline model up to degree 3.

The paper contributes to the factor zoo literature. cochrane2011presidential identifies the factor zoo phenomenon and calls for research to identify factors that are actually important for the SDF. harvey2016and argue that, based on a multiple-hypothesis-testing framework, many zoo factors are not significant. harvey2021lucky provide a bootstrap model selection framework that accounts for multiple hypothesis testing. Relatedly, bryzgalova2023bayesian and bryzgalova2025forest develop a Bayesian framework and a decision-tree framework for evaluating high-dimensional factor models, respectively. mclean2016does show that many zoo factors suffer from post-publication bias. green2017characteristics use Fama-MacBeth regressions to evaluate 94 firm characteristics and argue that fewer than ten provide significant independent information about expected returns. Our contribution to this literature is to show that a substantial fraction of the factor zoo can be explained by higher-order terms and interactions of common sparse linear factors.

Our paper is methodologically related to feng2020taming, but it addresses a different economic question. feng2020taming study a factor-screening problem: whether a proposed factor contributes incremental pricing information conditional on a large set of existing factors. We study, instead, whether much of the pricing information in the factor zoo can be accounted for by a small higher-order factor model built from a few common linear factors. Put differently, we ask whether the proliferation of factors largely reflects higher-order manifestations of a few common factors. Our forward selection Fama-MacBeth method is designed for this economic objective. Unlike the Lasso-based method in feng2020taming, which solves a penalized least-squares problem and induces sparsity through coefficient shrinkage, our procedure sequentially selects factors that maximize incremental cross-sectional pricing fit, thereby providing a natural ordered ranking of the importance of higher-order terms. This distinction is important both economically and algorithmically. First, forward selection does not shrink the risk-price estimates of selected factors, thereby preserving their economic interpretation. Second, the sequential nature of forward selection allows us to rank higher-order terms by their incremental contribution to pricing power. This ranking itself has economic content: for example, second-order size effects enter before cubic terms, and interactions dominate pure higher powers. Lasso does not provide such an ordered decomposition because variables can both enter and exit along the regularization path. Finally, our main theoretical contribution is to establish high-dimensional asymptotic theory for the forward selection Fama-MacBeth procedure. In particular, we derive the {\em rate of convergence} of the resulting estimator of SDF loadings in a setting where the number of candidate factors can be large relative to the time dimension. To the best of our knowledge, no such results are available for forward selection in the multi-step Fama-MacBeth framework. This theoretical contribution is conceptually distinct from feng2020taming, whose focus is on valid {\em inference} for testing the contribution of new factors using an $\ell_1$-regularized (Lasso-based) estimator.

Our paper also relates to a vast literature on factor models in asset pricing. A major part of the literature focuses on linear factor models, such as the characteristic-based models (e.g., fama1993common fama1993common,fama1996multifactor,fama2015five; Carhart1997, Carhart1997; koijen2017cross, koijen2017cross) and statistical-based models (e.g., connor1986performance, connor1986performance; sandulescu2021model; korsaye2023global). A small subset of papers emphasizes non-linearity in explaining asset returns, such as bansal1993no, harvey2000conditional, and dittmar2002nonlinear. We contribute to this literature by showing that higher-orders and higher-order interactions of common linear factors are important components of the SDF and price the cross-sectional asset returns.

Finally, our paper relates to the machine learning literature. gu2020empirical show that machine learning methods are better in predicting future stock returns out-of-sample than OLS methods. freyberger2020dissecting propose a group LASSO method to choose firm characteristics that are most related to future returns. kozak2020shrinking show that machine learning methods applied to the factor zoo can better approximate the SDF. lettau2020factors propose a method for estimating latent asset pricing factors that fit both the time-series and the cross-section of asset returns. giglio2021asset propose a method to estimate risk premia of factors even if there are missing factors. Borri2024factor show that a new non-linear single-factor asset pricing model motivated by the Kolmogorov-Arnold representation theorem is able to price assets from different asset classes. In this paper, we propose a forward selection method to estimate a high-dimensional stochastic discount factor model, isolating the most relevant higher-order factors and showing their importance.

The rest of the paper is organized as follows. Section (ref) discusses the forward selection Fama-MacBeth procedure. Section (ref) presents the data and the main empirical results. Section (ref) presents additional robustness analyses. The Online Appendix contains technical derivations and further empirical results.

Methodology

In this section, we set up the model. We introduce its main components -- returns, factors, and the stochastic discount factor -- and explain the relationship between them. We then present a forward selection Fama-MacBeth regression procedure as a method for selecting the factors that drive the stochastic discount factor. Using this selection procedure, we propose an estimator of the stochastic discount factor loadings. As we will see, the proposed estimator is biased. We therefore also explain how to debias it and how to perform inference on the stochastic discount factor loadings using the debiased estimator.

Notation. For any integer $p\geq 1$, we denote $[p]=\{1,\dots,p\}$ and $\mathbf 0_p = (0,\dots,0)^{\top}\in\mathbb{R}^p$. For any set $\mathrm{S}\subset[p]$, we define $\mathrm{S}^c = [p]\setminus\mathrm{S}$, and we denote by $|\mathrm{S}|$ the number of elements in $\mathrm{S}$. For any vector $x = (x_1,\dots,x_p)^{\top}$ and any set $\mathrm{S}\subset[p]$, we use $x_{\mathrm{S}} = (x_i)^{\top}_{i\in\mathrm{S}}$ to denote the sub-vector of $x$ consisting of all its components with indices in $\mathrm{S}$. Similarly, for any vector $x_t = (x_{t,1},\dots,x_{t,p})^{\top}$ with an extra index $t$ (or some other index), we denote by $x_{t,\mathrm{S}} = (x_{t,i})_{i\in\mathrm{S}}^{\top}$ the corresponding sub-vector. For any matrix $A = (A_{i,j})_{i\in[N],j\in[p]}$ and any set $\mathrm{S}\subset[p]$, we use $A_{\mathrm{S}} = (A_{i,j})_{i\in[N],j\in\mathrm{S}}$ to denote the sub-matrix of $A$ consisting of all its columns with indices in $\mathrm{S}$. For any vector $x = (x_1,\dots,x_p)^{\top}\in\mathbb{R}^p$, we let $\|x\|_0 = \sum_{j\in[p]} \mathbf I\{x_j\neq 0\}$ be its $\ell_0$-“norm” and let $\|x\|_1 = \sum_{j\in[p]} |x_j|$, $\|x\|_2 = (\sum_{j\in[p]} |x_j|^2)^{1/2}$, and $\|x\|_{\infty} = \max_{j\in[p]}|x_j|$ be its $\ell_1$-, $\ell_2$-, and $\ell_{\infty}$-norms. For any matrix $A = (A_{i,j})_{i\in[N],j\in[p]}\in\mathbb{R}^{N\times p}$, we let $\|A\|_2 = \sup_{v\in\mathbb{R}^p\colon \|v\|_2=1} \|Av\|_2$ and $\|A\|_{\infty,1} = \max_{i\in[N]}\sum_{j\in[p]} |A_{i,j}|$ be its spectral and $\ell_{\infty,1}-$norms. For any symmetric matrix $A\in\mathbb{R}^{p\times p}$, we let $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ be its smallest and largest eigenvalues. We use the abbreviation “\ensuremath{w.p.a.1}” to denote “with probability approaching one.” For non-random sequences $\{a_{T}\}_{T\geq1}$ and $\{b_{T}\}_{T\geq 1}$ in $\mathbb{R}$, we write $a_{T}\lesssim b_{T}$ if there exists a bounded non-random sequence $\{C_{T}\}_{T\geq 1}$ such that $a_{T}= C_{T} b_{T}$ for all $T\geq 1$. For random sequences $\{a_{T}\}_{T\geq 1}$ and $\{b_{T}\}_{T\geq 1}$ in $\mathbb{R}$, we write $a_{T} \lesssim_P b_{T}$ if there exists a bounded in probability random sequence $\{C_{T}\}_{T\geq1}$ such that $a_{T} = C_{T} b_{T}$ for all $T\geq 1$.

Model

Let $r_{i,t}\in\mathbb{R}$ be the excess return of asset $i\in[N]$ at time period $t\in[T]$ and let $f_t = (f_{t,1},\dots,f_{t,p})^{\top}\in\mathbb{R}^p$ be the vector of factors. Throughout the paper, we assume that the returns and the factors are stationary over time. In addition, we assume that the vector of factors is able to span the stochastic discount factor (SDF):

where $m_t$ is the SDF and $\psi = (\psi_1,\dots,\psi_p)^{\top}\in\mathbb{R}^p$ is the vector of SDF loadings.

By definition, since we work with the {\em excess} returns, $m_t$ is orthogonal to the returns in the sense that $\mathrm{E}[m_t r_{i, t}] = 0$ for all $i\in[N]$. Using this equation in combination with (ref), it then follows from Section 6.3 in C05 that for each asset $i\in[N]$,

where $\beta_{i} = (\beta_{i,1},\dots,\beta_{i,p})^\top\in\mathbb{R}^p$ is the vector of factor exposures corresponding to the factors $f_t$, $\gamma = (\gamma_1,\dots,\gamma_p)^{\top}\in\mathbb{R}^p$ is the vector of corresponding factor risk premia, and $\varepsilon_{i,t}\in\mathbb{R}$ is the residual satisfying $\mathrm{E}[\varepsilon_{i,t}] = 0$ and $\mathrm{E}[\varepsilon_{i,t}f_t] = \mathbf 0_p$. The vector of SDF loadings $\psi$ is related to the vector of factor risk premia $\gamma$ via the following equation:

where we assumed that the matrix $\mathrm{E}[(f_t - \mathrm{E}[f_t])(f_t - \mathrm{E}[f_t])^{\top}]$ is invertible. In addition, by using standard least-squares projection formulas together with (ref) and (ref), one can show that

where $\mathrm{C}_{i,j} = \mathrm{E}[(f_{t,j}-\mathrm{E}[f_{t,j}])(r_{i,t} - \mathrm{E}[r_{i,t}])]$ is the covariance between $r_{i,t}$ and $f_{t,j}$.

In this section, our task is to propose a new estimator of the vector of SDF loadings $\psi$ that can be used in the case when the number of factors $p$ is large, potentially much larger than the number of assets $N$ and/or the number of time periods $T$. As we will see, the classic estimation techniques may fail in this case. In fact, one can show that the SDF loadings cannot be consistently estimated without further assumptions in this case. Therefore, to make progress, and following the literature on estimation of high-dimensional models, we assume that the model is sparse, in the sense that most of the variation in the SDF in equation (ref) can be explained by a relatively small number of factors. In other words, only a relatively small number of coefficients in the vector $\psi$ are sufficiently different from zero.

Estimation

To develop intuition for the proposed method, let us first recall the classic Fama-MacBeth procedure that can be used to estimate the vector of risk premia $\gamma$ in the low-dimensional case, i.e. when $p$ is small. The procedure consists of two steps. First, using equation (ref), for each asset $i\in[N]$, we run a time series OLS regression of $r_{i,t}$ on $f_t$, with intercept included, to obtain estimates of factor exposures $\widehat\beta_i$. Second, we run a cross-sectional OLS regression of average return $\bar r_i = T^{-1}\sum_{t=1}^T r_{i,t}$ on $\widehat\beta_i$ (without intercept). The vector of slope estimates in this regression is then an estimator of the vector of risk premia $\gamma$. In turn, once the vector of risk premia is estimated, the vector of SDF loadings $\psi$ can be estimated using an empirical version of equation (ref).

The procedure described above is consistent when the number of factors $p$ is small. When the number of factors $p$ is large, however, all three steps may yield large estimation errors or fail as they require inverting estimates of poorly-conditioned/singular matrices. To deal with this problem, we propose a forward selection algorithm that starts with no factors and then adds factors iteratively one by one with the goal of achieving the highest possible $R^2$ in the cross-sectional regression of the Fama-MacBeth procedure. We stop the algorithm once a desired number of factors is achieved or the growth in the $R^2$ becomes sufficiently small. Once a set of factors is selected, we perform the Fama-MacBeth procedure and use an empirical version of equation (ref) to estimate $\psi$ as described above on selected factors. We set an estimate of the SDF loadings for non-selected factors to be zero.

To define our method more precisely, for each subset $\mathrm{S}$ of $[p] = \{1,\dots,p\}$, consider the Fama-MacBeth procedure using only factors in the set $\mathrm{S}$. In particular, define

where $\bar f = (\bar f_1,\dots,\bar f_p)^{\top} = T^{-1}\sum_{t=1}^T f_{t}$ and the inverse of the matrix is understood as the Moore-Penrose generalized inverse if the matrix is singular.\footnote{Note that the vector $\widehat\beta_{i,\mathrm{S}}$ is not a sub-vector of some $p$-dimensional vector $\widehat\beta_i$. In fact, even if we have $j\in\mathrm{S}_1\cap\mathrm{S}_2$ for some $j$, it is clearly possible that components of the vectors $\widehat\beta_{i,\mathrm{S}_1}$ and $\widehat\beta_{i,\mathrm{S}_2}$ corresponding to the factor $j$ may differ.} Also, let $R_{FM}^2(\mathrm{S})$ be the $R^2$ in the cross-sectional OLS regression of $\bar r_i$ on $\widehat\beta_{i,\mathrm{S}}$. We apply the forward selection algorithm with the aim to maximize $R_{FM}^2(\mathrm{S})$. For a desired number of factors $\widehat s$, the algorithm takes the following form:

Alternatively, the algorithm can be stopped when adding an extra factor does not significantly increase the $R^2$, i.e. $R_{FM}^2(\mathrm{S}\cup \widehat j) - R_{FM}^2(\mathrm{S}) \leq \epsilon$ for some small value of the tolerance parameter $\epsilon$ and $\widehat j$ appearing in Step 2 of the algorithm. In practice, we use $\epsilon = 1\%$, which seems meaningful from the empirical point of view.\footnote{Additionally, we implement a cross-validation method to select the optimal $\epsilon$ in Appendix (ref).} Also, some factors that are deemed to be important, e.g. Fama-French factors, can be included in the algorithm in the first step, i.e. we can set $\mathrm{S} = \mathrm{S}_0$ for some pre-specified set of factors $\mathrm{S}_0$ in Step 1 of the algorithm instead of setting $\mathrm{S} = \emptyset$.

The result of Algorithm (ref) is the set of factors $\mathrm{S}$, which we denote as $\widehat{\mathrm{S}}$ throughout the paper. Using this set, we estimate the vector of SDF loadings by $\widehat\psi = (\widehat\psi_1,\dots,\widehat\psi_p)^{\top}$, where we set

and $\widehat\psi_j = 0$ for all $j\in\widehat{\mathrm{S}}^c$. We will prove consistency and derive the rate of convergence of this estimator in Section (ref) below.

Debiased Estimation

Here, we explain how to debias the estimator $\widehat\psi = (\widehat\psi_1,\dots,\widehat\psi_p)^{\top}$ proposed in the previous section. To develop intuition for the debiasing procedure, fix any component of the vector $\widehat\psi$, say $\widehat \psi_j$, and observe that even though Algorithm (ref) is expected to select factors whose SDF loadings are sufficiently different from zero, it may not select factors whose SDF loadings are close to zero. The estimator $\widehat\psi_j$ may thus be biased for two reasons. First, $\psi_j$ may be close to zero itself, in which case Algorithm (ref) may not be able to select it, leading to the bias in $\widehat\psi_j$ as long as $\psi_j\neq 0$. Second, even if $\psi_j$ is not close to zero and the factor $j$ is selected, Algorithm (ref) may fail to select other factors $k$ whose covariance with returns $\mathrm{C}_{i,k}$ is cross-sectionally correlated with $\mathrm{C}_{i,j}$, leading to the omitted variable bias in $\widehat\psi_j$ as long as $\psi_k\neq 0$, which follows from the OLS characterization of the estimator $\widehat\psi_{\widehat{\mathrm{S}}}$ in equation (ref).

Two reasons for the bias of $\widehat\psi_j$ suggest that we can debias it by appropriately expanding the set $\widehat{\mathrm{S}}$. Specifically, we need to add to the set $\widehat{\mathrm{S}}$ the factor $j$ itself and also all factors $k$ whose covariance with returns $\mathrm{C}_{i,k}$ is cross-sectionally correlated with $\mathrm{C}_{i,j}$. In practice, since neither $\mathrm{C}_{i,j}$ nor $\mathrm{C}_{i,k}$ is observed, we work with sample covariances $\widehat \mathrm{C}_{i,j}$ and $\widehat \mathrm{C}_{i,k}$.

More formally, we proceed as follows. First, as before, let $\widehat{\mathrm{S}}$ be the set $\mathrm{S}$ produced by Algorithm (ref). Second, let $\widehat{\mathrm{S}}_{j}$ be the set produced by Algorithm (ref) with $R^2_{FM}(\mathrm{S})$ replaced by the $R^2$ in the cross-sectional OLS regression of $\widehat \mathrm{C}_{i,j}$ on $\widehat \mathrm{C}_{i,\mathrm{S}}$ (without intercept) and with $[p]$ replaced by $[p]\setminus\{j\}$. Third, set $\widehat{\mathrm{S}}_{D,j} = \widehat{\mathrm{S}}\cup \widehat {\mathrm{S}}_{j}\cup \{j\}$. Fourth, let $\widehat\psi$ be the estimator defined in Section (ref) but with the set $\widehat{\mathrm{S}}$ there replaced by the set $\widehat{\mathrm{S}}_{D,j}$ here. The corresponding component $\widehat\psi_j$ of the vector $\widehat\psi$ will be a debiased estimator of $\psi_j$, which we denote as $\widehat\psi_{D,j}$. By repeating this procedure for all $j\in[p]$, we obtain the full vector $\widehat\psi_D = (\widehat\psi_{D,1},\dots,\widehat\psi_{D,p})^{\top}$. In the next section, we will prove that for each $j\in[p]$, the estimator $\widehat\psi_{D,j}$ is $\sqrt T$-consistent and $\sqrt T(\widehat\psi_{D,j} - \psi_j)$ converges to a centered normal distribution.

Asymptotic Theory

Throughout the rest of the paper, we consider asymptotics where $T$ gets large, i.e. $T\to\infty$. We assume that the number of assets $N$ and the number of factors $p$, as well as all coefficients and distributions of random vectors in the model, are allowed to vary with $T$. In most cases, however, we keep this dependence implicit.

Let $k,K>0$ be constants such that $k\leq K$. Also, let $\bar s_{T}\geq 1$ be an integer that is allowed to depend on $T$ (and on $N$ and $p$ through $T$). In addition, let $\mathrm{C}=(\mathrm{C}_{i,j})_{i\in[N],j\in[p]}$ be the matrix of covariances between returns $r_{i,t}$ and factors $f_{t,j}$ and let $\widehat{\mathrm{C}}=(\widehat{\mathrm{C}}_{i,j})_{i\in[N],j\in[p]}$ be the corresponding matrix of sample covariances. Moreover, for brevity of notation, for all $t\in[T]$, let $v_t = (v_{t,1},\dots,v_{t,p})^\top = f_t - \mathrm{E}[f_t]$ be the demeaned version of the vector of factors $f_t$. Finally, to state some of the assumptions below in a compact form, let $u_{1,t} = 1$ and $u_{2,t} = \psi^{\top}v_t$ for all $t\in[T]$ and $u_{i+2,t} = \varepsilon_{i,t}$ and $u_{i+2+N,t} = \beta_{i}^{\top}v_t$ for all $i\in[N]$ and $t\in[T]$. Similarly, let $w_{t,1} = 1$ and $w_{t,2} = \psi^{\top}v_t$ for all $t\in[T]$ and $w_{t, j+2} = v_{t,j}$ for all $t\in[T]$ and $j\in[p]$.

Assumption (ref) is a quantitative version of the uniform law of large numbers. For example, given that $\mathrm{E}[\varepsilon_{i,t}v_{t,j}] = 0$, we expect from the classical central limit theorems for time series data that $T^{-1}\sum_{t=1}^T \varepsilon_{i, t}v_{t,j} \lesssim_P \sqrt{1/T}$. Provided that the random variables $\varepsilon_{i, t}v_{t,j}$ have sufficiently light tails, we then also expect that $T^{-1}\sum_{t\in[T]} \varepsilon_{i, t}v_{t,j} \lesssim_P \sqrt{\log (Np)/T}$ uniformly over $(i,j)\in[N]\times[p]$. Assumption (ref) imposes conditions of this type on various random variables appearing in our analysis below. Conditions like Assumption (ref) are very common in the literature on high-dimensional estimation.

Assumption (ref) states that the vector of SDF loadings $\psi$ is approximately sparse. In particular, it requires that a sub-vector $\psi_{\mathrm{S}_0}$ of the vector $\psi$ of size at most $\bar s_{T}$ is sufficient to explain most of the cross-sectional variation of mean returns; see equation (ref) for the relationship between SDF loadings and mean returns. Some researchers, e.g. nagel2021machine, argue against sparsity in finance but in different contexts, e.g. in the context of explaining cross-sectional variation of mean returns by past returns. We emphasize that the type of sparsity we impose is different. Indeed, {\em all} cross-sectional variation of mean returns is explained by just one factor: the stochastic discount factor. Thus, if it happens that one of the factors in the vector $f_t$ is equal to the stochastic discount factor, Assumption (ref) will be satisfied with a singleton set $\mathrm{S}_0$. Also, although we allow the set $\mathrm{S}_0$ to depend on $T$ (and on $N$ and $p$ through $T$), we omit this dependence for brevity of notation.

Provided $\lambda_{\min}(\mathrm{E}[v_tv_t^{\top}])$ is bounded away from zero and $\bar s_{T}$ does not grow too fast, Assumption (ref) can be proven to hold if the factors $f_t$ exhibit relatively weak time series dependence and have sufficiently light tails; e.g. see BC09.

The lower bound in Assumption (ref) is similar in flavor to that in Assumption (ref). It means that there is no multicollinearity in the vectors of covariances $\mathrm{C}_{\{j\}}$ when we restrict attention to a relatively small number of factors $j$. The upper bound should be viewed as a mild technical condition. Note also that Assumption (ref) is our key identification condition. Indeed, lewellen2010skeptical showed that without this assumption, it is possible that factors with zero SDF loadings can yield high $R^2$ in the cross-sectional OLS regression as long as they are correlated with the factors explaining the SDF. Assumption (ref) ensures that this phenomenon does not occur.

The first part of this assumption means that we carry out at most $\bar s_{T}$ rounds in our forward selection Fama-MacBeth procedure, Algorithm (ref), so that we select at most $\bar s_{T}$ factors. This condition is rather natural given the assumed sparsity of the true vector $\psi$. The second part of Assumption (ref) implies that the number of selected factors $|\widehat\mathrm{S}|$ should be larger than the number of important factors $|\mathrm{S}_0|$ at least by a multiplicative constant proportional to $\log T$. We impose this condition because the forward selection algorithm may occasionally select unimportant factors by chance. The third part of Assumption (ref) gives an upper bound on the number of selected factors but is less essential than the first two parts. Indeed, selecting more factors than the number in the upper bound will undermine the convergence rate of our estimator $\widehat\psi$ but will not make it inconsistent unless way too many factors are selected.

Assumption (ref) is a very mild regularity condition that is expected to hold in most cases. Assumption (ref) restricts how fast $\bar s_{T}$, $N$, and $p$ are allowed to grow relative to $T$. Most importantly, it allows the number of factors $p$ to be much larger than the number of time periods $T$. Together with Assumption (ref), it also requires that the size of the approximating vector $\psi_{\mathrm{S}_0}$ is small relative to the number of time periods $T$.

To formulate the asymptotic normality result for the debiased estimator $\widehat\psi_D$, for all $j\in[p]$, let $\eta_j = (\eta_{j,1},\dots,\eta_{j,p})^{\top}\in\mathbb{R}^p$ be the vector defined by $\eta_{j,j} = 0$ and

so that $\eta_{j,\{j\}^c}$ is the vector of coefficients of the time-series least-squares projection of the demeaned factor $v_{t,j}$ on other demeaned factors $v_{t,\{j\}^c}$. Also, for all $t\in[T]$, let $ z_{t,j} = v_{t,j} - \eta_{j,\{j\}^c}^{\top}v_{t,\{j\}^c} $ be the corresponding residual and let $\sigma_{z,j}^2 = \mathrm{E}[z_{t,j}^2]$ be the residual variance. The asymptotic normality result requires several additional assumptions but, for brevity of the main text, we have placed them in Appendix (ref).

Empirical Analysis

In this section, we present the data used in the empirical analysis and the main results of the selection of the higher-order factors by the forward selection Fama-MacBeth procedure.

Data

We present a description of the data used for the estimation of the FS-FMB procedure. In order to estimate the model and perform cross-sectional asset pricing tests, we use data from various sources. Test assets are 484 characteristic-managed portfolios from kozak2020shrinking, ranging from October 1973 to December 2019.\footnote{ Portfolios are from Serhiy Kozak's web page (\href{https://serhiykozak.com/\#/data-code}{https://serhiykozak.com/\#/data-code}). } These portfolios are constructed by weighting stocks based on their characteristic signals, which are based on cross-sectional ranks. Small firms (below 0.01% of market cap) are excluded. Baseline factors are the fama2015five's five factors and the momentum factor of \citet*{jegadeesh1993returns}, which we obtain from Kenneth French's data library. We refer to these six baseline factors as the FF5M model. Furthermore, we consider factors for 148 characteristics in 13 themes, using data from the US equity market, constructed by \citet*{jensen2023there}.\footnote{Factors are capped value weighted from Global Factor Data (\href{https://jkpfactors.com/}{https://jkpfactors.com/}).} We refer to these 148 factors as the factor zoo (see, e.g., \citet*{feng2020taming}).

We construct factors based on the FF5M factors. We denote both the higher-orders and the interactions of FF5M as higher-order factors. Specifically, for any factor $f_{i}$ in the FF5M, we first form all the $f_{i}^{2}$ and $f_{i}^{3}$ as candidate factors. Moreover, we form all the pairwise interactions of degree 2 and 3. That is, for any two factors $f_{i}$ and $f_{j}$ in the FF5M model, we form the pairwise interactions of the type $f_{i}\times f_{j}$ and $f_{i}\times f_{j}^{2}$ as candidate factors. Therefore, the total number of candidate higher-order factors is 57, where there are 6 $f_{i}^{2}$, 6 $f_{i}^{3}$, 15 $f_{i}\times f_{j}$, and 30 $f_{i}\times f_{j}^{2}$.

Table (ref) presents summary statistics of the test assets and factors. For the mean, the standard deviation, and the quantiles, we report averages across assets. All returns are monthly, in excess of the US risk-free rate, and reported in percentage. The average mean return for the test assets is 0.65% monthly, with a standard deviation of 4.97%. The average return for the FF5M factors is lower and equal to 0.39% monthly, with a standard deviation of 3.16%. The average return for the factor zoo is even lower and equal to 0.23% monthly, with a standard deviation of 2.71%. The higher-order factors have a distribution characterized by more extreme realizations, as exemplified by the 95% quantile of the powers of degree 2 and 3 equal, respectively, to 42.75% and 167.67%. We note that the higher-order factors are generally non-tradable. Figure (ref) in the Appendix plots the distributions of the returns of the higher-order factors.

Selected Factors by Forward Selection Fama-MacBeth Regressions

Our goal is to study whether and which of the higher-order factors constructed from the six factors in the baseline FF5M are important in explaining the SDF. To select the most relevant higher-order factors to the SDF, we use the FS-FMB procedure that we have discussed in detail in Section (ref) above. We start with the six baseline FF5M factors and add the higher-order interactions one at a time.

We present the main estimation results from the second step of the FS-FMB procedure in Table (ref). Panel A shows the performance of the baseline models using the 484 characteristic-managed portfolios from kozak2020shrinking. The baseline models include the standard CAPM, the fama1993common's 3-factor model (FF3), the fama2015five's 5-factor model (FF5) and the 5-factor model augmented with the momentum factor of jegadeesh1993returns (FF5M). The adjusted cross-sectional R-squared for the CAPM is just 0.032, and the estimated intercept (0.009) is significantly different from zero with a $t$-statistic of 3.6. The adjusted cross-sectional R-squared reaches 0.312 for the FF5M model, while the estimate for the intercept (0.004) remains significantly different from zero with a $t$-statistic of 2.2.

Panel B refers to the FS-FMB procedure. At each step, we augment the FF5M by one of the additional higher-order factors, which we select in order to maximize the adjusted cross-sectional R-squared. The selection procedure stops when the improvement in the R-squared in a given step is smaller than 1 pp. The forward selection FMB procedure selects 7 higher-order factors: the SMB factor squared (SMB2), the interaction between SMB2 and the momentum factor (SMB2{*}Mom), the interaction between the momentum factor squared and the profitability factor (Mom2{*}RMW), the market factor squared (Mkt-RF2), the interaction between Mkt-RF2 and the profitability factor (Mkt-RF2{*}RMW), the interaction between the market and size factors (Mkt-RF{*}SMB) and the interaction between the value factor squared and the market factor (HML2{*}Mkt-RF). Of these 7 higher-order factors that are selected, 2 are higher-orders (SMB2 and Mkt-RF2) and 5 are interactions (SMB2{*}Mom, Mom2{*}RMW, Mkt-RF2{*}RMW, Mkt-RF{*}SMB, and HML2{*}Mkt-RF). Interestingly, none of the higher-orders of degree 3 are selected.

The adjusted cross-sectional R-squared increases from 31.2% for the baseline FF5M model to 41% in the first step when SMB2 is included, and gradually increases to 59%, in the last step when all 7 factors are included and the procedure stops. Across the seven steps, the estimate for the intercept is around 0.003, with a $t$-statistic in the final step of 1.819, which is below the threshold for significance at the 5% level.

We then examine the SDF loadings of the higher-order factors selected by the FS-FMB procedure. The results are summarized in Table (ref). Panel A shows the non-debiased estimates based on the standard FMB two-pass regression, while Panel B presents the debiased estimates as outlined above. The directions of the point estimates based on the non-debiased procedure and the debiased procedure are the same. Moreover, the estimates from the two procedures are similar in magnitude. The estimates based on the non-debiased procedure are all significant at the 5% level, while the estimates based on the debiased procedure are all significant at the 10% level with 5 of the 7 SDF loading estimates being significant at the 5% level. Overall, the results show that the higher-order factors are important components of the SDF and significantly price assets in the cross-section.

Figure (ref) plots predicted against realized returns for all the test assets. Panel A shows the performance of the FF5M model and panel B presents the higher-order model. The FF5M has a limited ability in explaining the average returns (Panel A), in line with the relatively low R-squared of the model. On the other hand, there is a strong positive relationship between the average realized returns and the predicted returns for the higher-order model selected by the FS-FMB procedure.

We note that the predicted return used in the graph is the OLS estimate of the betas from the first step of the two-pass FMB procedure times the estimated prices of risk. That is, we do not impose the no-arbitrage condition for the tradable factors that pins down the price of risk by the sample mean of the factor. Doing so would, by construction, result in an even worse fit for the baseline FF5M models where all the factors are tradable. We present the results in Table (ref) in the Online Appendix where we impose the restriction for the tradable factors. The results confirm the importance of higher-order factors in pricing assets in the cross-section.

Out-of-Sample

The main results presented in Table (ref) are based on in-sample estimation. We further provide evidence of the out-of-sample performance of the higher-order model selected by the FS-FMB procedure. We consider two types of out-of-sample tests: (1) out-of-sample in terms of asset space and (2) out-of-sample in the time dimension.

Asset space

In Table (ref), we present an additional set of results based on cross-validation applied to the FS-FMB procedure. The cross-validation is applied to the second step of the FS-FMB procedure. This is an out-of-sample exercise in terms of asset space, where we select higher-order factors and estimate SDF loadings from the training assets, and then evaluate the out-of-sample performance in the out-of-sample assets.

Specifically, at each step of the FS-FMB procedure, we divide the sample into 5 equal-sized, non-overlapping folds. For each fold $k$, we use the remaining 4 folds as the training set and the left-out fold as the test set. In each training step, we estimate the cross-sectional regression of average returns on the betas and obtain estimates for the risk prices. Next, we use the estimated risk prices to compute the adjusted R-squared on the test set, where the degrees-of-freedom adjustment uses the number of test assets in each fold. This process is repeated 5 times and the adjusted cross-sectional R-squared associated with a given higher-order factor is the average, at each step, of the adjusted cross-sectional R-squareds obtained in the 5-fold cross-validation. The procedure stops when the improvement in the cross-validated adjusted R-squared falls below 1 pp. Table (ref) reveals that the FS-FMB procedure with cross-validation (CV) selects a total of 7 higher-order factors, with an in-sample adjusted cross-sectional R-squared in the final step of 0.582. The selected factors are the same as the first 6 from the baseline method. The cross-validated adjusted cross-sectional R-squareds, which represent out-of-sample measures of fit of the cross-sectional model, reach 0.471 in the final step. Overall, these results show that the selected higher-order factor model has stable out-of-sample performance. In Appendix (ref), we further show that the results are robust to the choice of stopping threshold by selecting it optimally and entirely out-of-sample via nested cross-validation (Table (ref)).

Time periods

We next evaluate the out-of-sample performance in the time-series of the higher-order factor model against the benchmark factor models (i.e., CAPM, FF3, FF5, and FF5M). Specifically, we first split the sample into two 50-50 subsamples based on the time series, where we denote the first half of the sample as the training sample and the second half as the out-of-sample. Then, we estimate the covariances and the SDF loadings in the training sample based on the different factor models. We denote the $R^{2}$ from these estimations as $R_{train}^{2}$. Lastly, we compare the average asset returns in the out-of-sample to the predicted asset returns using covariances and the SDF loadings estimated from the training sample to obtain the out-of-sample $R^{2}$, or $R_{oos}^{2}$. We note that both the $R_{train}^{2}$ and the $R_{oos}^{2}$ are likely to be lower than the full sample $R^{2}$ for any model due to estimation errors, especially related to average returns. It is well-known that estimations of the average equity returns are noisy and require a long time series. Therefore, our focus is comparing $R_{oos}^{2}$ for the higher-order factor model against those of the benchmark factor models.

Table (ref) presents the results. The $R_{oos}^{2}$s of CAPM and the FF3 factor model are close to zero. The $R_{oos}^{2}$s of the FF5 factor model and the FF5M factor model are 7% and 8.9%, representing an improvement from the CAPM and the FF3 factor model. The $R_{oos}^{2}$ of the higher-order factor model selected from the FS-FMB procedure is significantly larger than the benchmark factor models at 14.8%, which is about twice as high as the $R_{oos}^{2}$s of FF5 and FF5M factor models. The $R_{oos}^{2}$ of the higher-order factor model where higher-order factors are selected in the training sample is similar at 15.8%.

One concern is that the 50-50 split where the first half of the sample is the training sample and the second half of the sample is the out-of-sample may be arbitrary. As a robustness check, we repeat the analyses by randomly selecting, for 1,000 times, 50% of the sample as the training sample and the remaining as the out-of-sample. We report the results in Table (ref) in the Online Appendix. Table (ref) delivers the same message -- the $R_{oos}^{2}$ of the higher-order factor model is significantly higher than those of the benchmark factor models.

Factor Zoo

In this subsection, we show that the baseline FF5M model, augmented with the higher-order factors selected by the FS-FMB procedure, substantially reduces the majority of the factors in the factor zoo, or zoo factors.

We employ two empirical methods to support this argument. First, we perform FMB regressions to test the significance of the SDF loading of each factor in the factor zoo, controlling for our higher-order factors. Second, we construct factor mimicking portfolios for the seven higher-order factors to conduct time-series asset pricing tests on the zoo factors.

Cross-sectional method

The first empirical strategy to evaluate whether the higher-order factors can reduce the zoo factors is based on the two-pass Fama-MacBeth procedure, motivated by the analysis in feng2020taming. Specifically, for each zoo factor, we estimate its SDF loading while controlling for our higher-order factor model, which includes the FF5M and the selected higher-order factors.

Figure (ref) summarizes the results. Panel A plots the distribution of the absolute values of the t-statistics for the intercept from the cross-sectional regressions. Panel B displays the distribution of the absolute $t$-statistics associated with the SDF loading of the factor zoo, using newey1987simple corrected standard errors. In both panels, the blue histogram represents the estimates obtained by controlling for the FF5M factors, while the green histogram represents the estimates obtained by additionally controlling for the higher-order factors selected by the forward selection FMB procedure. Vertical dashed lines indicate the median values of each distribution -- a vertical blue dashed-line for estimates obtained controlling for the FF5M factors and a vertical green dashed-line for estimates obtained controlling for the higher-order factor model. Additionally, we also include a vertical red dashed-line which corresponds to the median value with no controls, or the factor model including only the factor in the factor zoo.

Panel A shows that the intercept estimates tend to be highly statistically significant, with a median value of 4.00, when the factor model includes only the factor in the factor zoo. When the factor model includes both FF5M and the factor in the factor zoo, the median $t$-statistic decreases substantially but remains significant at the 5% level at 2.23. When the factor model further includes the higher-order factors, the median $t$-statistic is no longer statistically significant at the 5% level at 1.84.

Panel B shows that, when the factor model includes only the factor in the factor zoo, the median SDF loading estimate for the factor zoo is not statistically different from zero, with a $t$-statistic of 0.86. Note that these estimates are likely biased because of omitted factors (see, e.g., feng2020taming). Controlling for the FF5M, the median $t$-statistic increases substantially to 1.38. 34% of zoo factors become statistically significant at the 5% level with a $t$-statistic greater than 2 suggesting that FF5M helps to remove some bias in the factor estimates. Importantly, when we additionally control for the higher-order factors, we find that for 95% of the zoo factors the risk price estimate is not statistically different from zero, with the median $t$-statistic dropping from 1.38 to only 0.71. Specifically, we find that only 7 out of 148 zoo factors remain significant at the 5% level. Of these seven factors, three belong to the “low leverage” category of jensen2023there (cash_at, rd5_at, and rd_sale). The remaining significant factors are in the categories “size” (rd_me), “seasonality” (seas_6_10an), “accruals” (cowc_grla) and “quality” (\textit{qmj}). Overall, these results reveal that the higher-order factors selected by the forward selection FMB procedure substantially reduce the majority of the zoo factors.

Time-series method

The second empirical strategy we employ to evaluate whether the higher-order factors can capture the factor zoo is based on time-series regressions, where we study whether the zoo factors can be spanned by the FF5M and the selected higher-order factors.

Note that although the FF5M factors are tradable, the higher-order factors are generally not directly tradable. Therefore, in order to estimate time-series spanning regressions, we must first construct factor mimicking portfolios of the higher-order factors based on tradable assets. In order to do so, we project each higher-order factor onto the entire factor zoo and obtain its factor-mimicking portfolio as the regression's predicted values. These factor-mimicking portfolios only capture the information in the factor zoo that is related to the higher-order factors. The benefit of the factor-mimicking portfolio approach is that we are able to convert the non-tradable higher-order factors to tradable factors, but we note that the factor-mimicking portfolios of the higher-order factors inevitably lose information from the original factors. Figure (ref) shows the adjusted $R^{2}$ of projecting each selected higher-order factor onto the factor zoo. Mkt-RF{*}SMB has the lowest adjusted $R^{2}$ which is less than 40% and SMB2{*}Mom has the highest adjusted $R^{2}$ at more than 80%, suggesting the projections capture important fractions, though to various degrees, of variations in the higher-order factors. Nevertheless, we use them to check if they capture pricing information of each factor.

To do so, we estimate a time-series regression for each factor in the factor zoo on the FF5M and the factor-mimicking portfolios of the higher-order factors. The object of interest is the intercept. If the higher-order factors indeed span a given factor in the factor zoo, we expect to find a small and statistically insignificant intercept.

Figure (ref) summarizes the results of the time-series spanning regressions. Panel A plots the distribution of the absolute $t$-statistics associated with the intercepts, and Panel B shows the distribution of the absolute magnitudes for the intercepts. The latter are annualized and expressed in percentage points.

Panel A shows that the median of the absolute $t$-statistic of the average returns for the zoo factors is statistically significant at the 5% level at 2.06. Controlling for the FF5M brings the median absolute $t$-statistic from 2.06 to 1.87, a decline of about 9%. Controlling for the higher-order factors brings the median absolute $t$-statistic further down to 1.49, a decline of about 21% from 1.87. Panel B reveals a similar pattern. The median of the absolute intercept estimate of the zoo factors drops from 2.92 pp to 1.69 pp controlling for the FF5M model. It further decreases to 1.34 pp when the higher-order factors are included.

Relationship between factor zoo and higher-order factors

In the previous subsections, we show that the higher-order factor model selected by the FS-FMB procedure is successful in accounting for a substantial amount of the zoo factors. Mechanically, the ability of the selected higher-order factor model to capture the factor zoo has to come from correlation between the zoo factors and some or all of the factors in the higher-order factor model. In this subsection, we provide information on the loadings of the zoo factors on the factors of the selected higher-order factor model, which sheds light on the relative importance of the higher-order factors in accounting for the zoo factors.

Table (ref) documents the fraction of zoo factors with loading significantly different from zero at the 5% confidence level on each of the FF5M factors and the higher-order factors selected by the forward selection FMB procedure. Column (1) shows the fraction of significance for the factors in FF5M model only and Column (2) reports the fraction of significance for the factors in the full higher-order factor model. Column (1) shows that the fractions of significance range from 45.9% for Mkt-RF to 75% for RMW for the 6 factors in FF5M. Column (2) shows that the fractions of significance remain stable for the 6 factors in FF5M, when we further include the selected higher-order factors. The fractions of significance for the higher-order factors range from 27% for Mkt-RF2 to 55.4% for HML2{*}Mkt-RF, indicating that it is not the case that a specific higher-order factor accounts for all the pricing results.

Higher-Order Factors and Macro Risk

In this subsection, we conduct analyses to better understand the economic forces the higher-order factors capture. In the first exercise, we examine the correlations between the selected higher-order factors and two important macroeconomic risks as proxied by the U.S. equity market returns and the CBOE VIX index.

In Table (ref), we present the results. Panel A documents the correlations with U.S. equity market returns and Panel B shows the correlations with CBOE VIX index. Beyond correlations over the full sample, we also calculate the correlations during the NBER recession periods and the correlations corresponding to the bottom 10%, the middle 10%-90%, and the top 10% of higher-order factors' return distribution.

For the first two higher-order factors (i.e., SMB2 and SMB2{*}Mom), the magnitude of the correlations with U.S. equity market returns is markedly higher during the NBER recession periods relative to the full sample. For example, the correlation between SMB2 and U.S. equity market returns is only 0.016 for the full sample, but it jumps to 0.205 during the NBER recession periods. For the last four higher-order factors (i.e., Mkt-RF2, Mkt-RF2{*}RMW, Mkt-RF{*}SMB, and HML2{*}Mkt-RF), they tend to be highly correlated with the VIX index, especially during the NBER recession periods. For example, the correlation between Mkt-RF2{*}RMW and VIX is 0.32 for the full sample and increases to 0.515 during the NBER recession periods. These results show that the selected higher-order factors capture different macroeconomic risks, especially during economic downturns.

The results in Table (ref) indicate that the higher-order terms---the powers and interactions---may systematically capture the state-dependent macroeconomic tail risks. To formalize this idea, in Panel A of Table (ref), we first evaluate the unconditional exposures of the selected higher-order factors to a broad set of common macroeconomic variables using OLS regressions. The results show that these macroeconomic variables are able to account for a meaningful fraction of the variation in the selected higher-order factors. The adjusted $R^{2}$ values range from 3% to 45%, exceeding 20% for five of the seven higher-order factors. Among these macroeconomic variables, measures related to aggregate uncertainty, in particular financial uncertainty, are highly associated with these higher-order factors.

While the baseline unconditional regressions establish a link to macroeconomic risk, if the higher-order terms capture tail dependence, we would expect the relationship to be more pronounced in the tails of the distribution. To better capture the potential state dependence, we estimate quantile regressions of the higher-order factors on the same set of macroeconomic variables. In Panel B of Table (ref), we report the conditional loadings for the bottom decile ($\tau=0.10$) and the top decile ($\tau=0.90$) of the higher-order factors on the macroeconomic variables. The extreme state dependence of these higher-order factors is most evident from the pseudo-$R^2$ of these quantile regressions. The explanatory power based on the pseudo-$R^2$ values for the tail regressions in Panel B of Table (ref) frequently exceeds 50%. Moreover, the coefficient estimates on financial uncertainty and measures of intermediary capital constraints become significantly more pronounced during these tail realizations.

These results suggest that the higher-order factors may capture non-linear exposure to macroeconomic tail risks, particularly those associated with financial intermediary uncertainty and distress. This interpretation is consistent with intermediary asset pricing models (e.g., he2013intermediary,brunnermeier2014macroeconomic), in which asset prices are disproportionately sensitive to intermediary stress during periods of extreme macroeconomic environments. The fact that the relationship between higher-order factors and macroeconomic variables, in particular financial intermediary constraints, strengthens precisely in extreme macroeconomic environments provides an economic rationale for why these non-linear terms of common factors are priced and why they may account for the broader factor zoo.

Additional Results

In this section, we provide additional results for our empirical analyses, and we discuss extensions of our baseline setup.

Simulation Results

The selected higher-order model incorporates 7 higher-order factors beyond the traditional FF5M model. Our results indicate that the model performance significantly improves with the 7 higher-order factors where the R-squared of the second-pass regression increases from 31.2% to 58.7%. However, as noted by lewellen2010skeptical, spurious factors can artificially inflate model fit in the second-pass of the Fama-MacBeth procedure, particularly when the test assets exhibit a strong factor structure. This concern is closely related to the weak-factor identification problem in linear asset-pricing models bryzgalova2015spurious. lewellen2010skeptical illustrate this concern using the Fama-French 25 Size-BM portfolios which have a particularly strong factor structure. While our analysis mitigates this issue by employing a broad set of 484 test assets and demonstrating the superior out-of-sample performance, we further conduct simulation exercises to confirm that our results are not driven by spurious factors.

Therefore, we estimate the sampling distribution of cross-sectional R-squareds obtained using randomly generated factors, following the approach advocated by lewellen2010skeptical. Specifically, in each simulation, we generate a set of 57 random factors, matching the number of candidate higher-order factors. These random factors are assumed to be $i.i.d.\sim N(0,\sigma^{2})$, and we set $\sigma^{2}$ equal to the sample variance of the market factor (Mkt-RF). We then apply a forward selection FMB procedure, sequentially selecting factors that maximize the adjusted cross-sectional R-squared when added to the FF5M factors. We repeat the procedure 1,000 times by generating 1,000 independent sets of factors.

The green histogram in Figure (ref) summarizes the sampling distribution of the resulting R-squareds from the 1,000 simulations. In the simulation exercise, we obtain a cross-sectional R-squared larger than 0.59 only in 0.1% of the simulations, meaning that the adjusted R-squared from the higher-order model is significant at the 0.1% level. This simulation exercise highlights the fact that it is extremely unlikely that our results are obtained by chance.

Note that the procedure can select a number of factors different than 7, which is the number of higher-order factors presented in Table (ref). In fact, in the simulation exercise we conduct, the maximum number of random factors selected is 12. In Panel A, we further compare the sampling distribution of R-squareds with that obtained by directly adding 7 randomly generated factors, $i.i.d.\sim N(0,\sigma^{2})$, to the FF5M factors, repeating the procedure 1,000 times (blue histogram). In this case, the maximum adjusted cross-sectional R-squared is only 0.463.

Finally, we compare the sampling distribution of R-squareds with that obtained with the additional constraint that the number of selected factors cannot be larger than 7 (blue histogram) in Panel B. In this case, we never observe R-squareds larger than 0.59 (maximum adjusted cross-sectional R-squared of 0.565).

An alternative factor zoo

In order to investigate the robustness of our results in terms of accounting for the factor zoo, we replicate the analysis with the two-pass Fama-MacBeth procedure using an alternative set of zoo factors, based on the factors constructed by ChenZimmermann2021. Specifically, for each factor in ChenZimmermann2021, we estimate its SDF loading while controlling for our higher-order factor model, which includes the FF5M and the selected higher-order factors. Figure (ref) summarizes the results. Panel A shows that the intercept estimates tend to be highly statistically significant, with a median value of 3.97, when the factor model includes only the factor in the factor zoo. When the factor model includes both FF5M and the factor in the factor zoo, the median $t$-statistic decreases substantially but remains significant at the 5% level at 2.23. When the factor model further includes the higher-order factors, the median $t$-statistic is no longer statistically significant at the 5% level at 1.87. Panel B reveals that, when the factor model includes only the factor in the factor zoo, the median SDF loading estimate for the factor zoo is not statistically different from zero, with a $t$-statistic of 0.80. Note that these estimates are likely biased because of omitted factors (see, e.g., feng2020taming). Controlling for the FF5M, the median $t$-statistic increases substantially to 1.91. 46% of zoo factors become statistically significant at the 5% level with a $t$-statistic greater than 2 suggesting that FF5M helps to remove some bias in the factor estimates. Importantly, when we additionally control for the higher-order factors, we find that for 93% of the zoo factors the risk price estimate is not statistically different from zero, with the median $t$-statistic dropping from 1.91 to only 0.72. Specifically, we find that only 12 out of 159 zoo factors remain significant at the 5% level. These results confirm our findings based on the factors in jensen2023there.\footnote{The factors in ChenZimmermann2021 that remain significant at the 5% confidence level are: MomSeason06YrPlus, CBOperProf, ChInv, STreversal, DelFINL, DivYieldST, InvGrowth, GrSaleToGrInv, MomOffSeason16YrPlus, DebtIssuance, CompositeDebtIssuance, OrderBacklog.}

Table (ref) in the Online Appendix presents the fraction of zoo factors in ChenZimmermann2021 with loading significantly different from zero at the 5% confidence level on each of the FF5M factors and the higher-order factors selected by the forward selection FMB procedure. Column (1) shows the fraction of significance for the factors in FF5M model only and Column (2) reports the fraction of significance for the factors in the full higher-order factor model. Column (1) shows that the fractions of significance range from 42.8% for CMA to 57.2% for SMB and RMW for the 6 factors in FF5M. Column (2) shows that the fractions of significance remain stable for the 6 factors in FF5M, when we further include the selected higher-order factors. The fractions of significance for the higher-order factors range from 20.8% for Mkt-RF2{*}RMW to 47.8% for SMB2{*}Mom. These results are consistent with the findings we obtain using the zoo factors in jensen2023there, further highlighting that it is not the case that a specific higher-order factor accounts for all the results.

Alternative sets of candidate factors

We start by considering alternative sets of higher-order factors formed using the FF5M factors. Table (ref) presents the estimation results of the forward selection FMB procedure with four alternative sets of higher-order factors. Specifically, Panel A presents estimation results using a set of higher-order factors of degree 2. The latter is a lower degree than the one from the baseline analysis (i.e., degree 3). These higher-order factors include all pairwise interactions of the FF5M, and the FF5M factors squared. In this case, the FS-FMB procedure selects four higher-order factors, reaching an adjusted cross-sectional R-squared of 0.517 in the final step, approximately 7 pp smaller than in the baseline setup (see Table (ref)). Panel B presents estimation results using a set of higher-order factors of higher degree, and specifically of degree up to 4. These factors, augment the higher-order factors of degree 3 of the baseline setup with the FF5M to the fourth power (i.e., $f_{i}^{4}$), and pairwise interactions of the FF5M factors of degree 4, i.e., $f_{i}^{2}\times f_{j}^{2},f_{i}^{3}\times f_{j},f_{i}\times f_{j}^{3}$. In this case, the FS-FMB procedure selects 9 higher-order factors, reaching an adjusted cross-sectional R-squared of 0.6 in the final step. The latter is a R-squared value very similar to that obtained in the baseline setup, with higher-order factors up to degree 3. These results suggest that higher-order factors up to degree 3 are most relevant in explaining the cross-section of asset returns, corroborating our choice of using higher-order factors of degree 3 in the baseline analysis.

Furthermore, Panel C and D report the results when augmenting the FF5M with only pairwise interactions of degree 2 and 3 (Panel C) and only the powers of degree 2 and 3 (Panel D). The results from these two panels reveal that incorporating pairwise interactions between FF5M is crucial. In fact, the model with interactions of degree 2 and 3 reaches an adjusted cross-sectional R-squared of 0.542, selecting 7 higher-order factors. In contrast, the model with powers of degree 2 and 3 selects 3 higher-order factors, for an adjusted cross-sectional R-squared of 0.467 in the final step. These results highlight that both the higher-order terms and the interactions of FF5M are important components of the SDF, especially higher-order interactions.