A penalized two-pass regression to predict stock returns with time-varying risk premia

A penalized two-pass regression to predict stock returns with time-varying risk premia

Keywords: two-pass regression, predictive modeling, large panel, factor model, LASSO penalization.\\

JEL classification: C13, C23, C51, C52, C53, C55, C58, G12, G17.

{\scriptsize$^{a}$Emlyon Business School, $^{b}$Geneva School of Economics and Management, University of Geneva, $^{c}$Faculty of Science, University of Geneva, $^{d}$Swiss Finance Institute.}{\scriptsize }

Introduction

Under the arbitrage pricing theory ross2013arbitrage, Chamberlain_Rothschild_1983, we know that risk premia are drivers of expected excess returns. Hence, estimating them should be useful for prediction of future equity excess returns. The workhorse to estimate equity risk premia in a linear multi-factor setting is the two-pass cross-sectional regression method developed by jensen1972capital and fama1973risk. A series of papers address its large and finite sample properties for linear factor models with time-invariant coefficients; see, for example, shanken1985multivariate,shanken1992estimation, jagannathan1998asymptotic, shanken2007estimating, kan2013pricing, and the review paper of jagannathan2010analysis (see bryzgalova2019bayesian for a recent Bayesian approach). In a time-varying setting, gagliardini2016time (henceforth referred as \citetalias{gagliardini2016time}) study how we can infer the dynamics of equity risk premia from large stock return data sets under conditional linear factor models (see also gagliardini2019estimation for a review of estimation of large dimensional conditional factor models in finance). They show how to explicitly account for the no-arbitrage restrictions relating the time-varying intercept and the time-varying factor loadings when writing the underlying linear regression to be estimated. In conditional factor models, we quickly loose parsimony in terms of covariates because of the cross-products induced by the no-arbitrage restrictions. chaieb2020factors show that a direct application of the \citetalias{gagliardini2016time} methodology in an international setting is challenging because of the large number of parameters needed to model the time-variations in factor exposures and risk premia. Applying the \citetalias{gagliardini2016time} methodology off-the-shelf to an international setting results in few or even zero stocks kept for several countries. To address this issue, they suggest to rely on iteratively selecting for each stock the most important covariates driving the dynamics of the factor loadings without violating the no-arbitrage restrictions.

The aim of this paper is to tackle this issue via LASSO-type penalisation techniques tibshirani1996regression to enforce sparsity for the time-variation drivers while also maintaining compatibility with the no-arbitrage restrictions. The shrinkage targets the time-invariant counterpart of the time-varying models. In a conditional factor setting, we aim at addressing the “multidimensional challenge” of cochrane2011presidential, namely select characteristics which really provide independent information about average excess returns. More specifically, the penalized first-pass (time-series) regression selects and estimates the regression coefficients ensuring a model specification compatible with the no-arbitrage restrictions through the Overlap Group-LASSO (OGL) of jacob2009group and its adaptive version of percival2012theoretical, the aOGL, which extends the original Group-LASSO of yuan2006model to groups of variables that may overlap. Indeed, if we do not introduce a quadratic term (or cross-products) in the time-varying intercept while the covariate is present in the time-varying factor loadings, we introduce ex-ante a model with arbitrage gagliardini2019estimation. By definition, we cannot estimate a coefficient for which its covariate is absent. On the contrary, if we delete a covariate in the time-varying factor loadings and keep it in the time-varying intercept, then its corresponding coefficients could be shrunk to zero by a standard LASSO for the first-pass regression, and thus could avoid ex-post a model with arbitrage if the true model is sparse. In a standard Ordinary Least Squares (OLS) first-pass procedure, those time-varying intercept coefficients could be estimated close to zero if the true model does not include that covariate in the time-varying factor loadings. By introducing groups based on finance theory derived from assuming no asymptotic arbitrage opportunities in the economy, our aOGL approach can only consider models compatible ex-ante with the no-arbitrage restrictions by construction. The groups take explicitly into account the links between the time-varying intercept and the time-varying loadings induced by the no-arbitrage restrictions. With only models satisfying ex-ante the no-arbitrage restrictions, we can substantially reduce the set of possible models within our model selection procedure. We derive an upper bound, and show that the number of possible models without grouping is divided by $2^3$, at least, and often by a much larger number in empirical applications. As an example, for the model specifications with four factors used in Section (ref), the set of possible models satisfying ex-ante the no-arbitrage restrictions is $2^{97}$ times smaller than the set of possible models without grouping. We exemplify this reduction with a simple two-factor example in Section (ref). It echoes the discussion in giannone2021economic that, if a prediction model with many predictors “lacks any additional structure, then there is no hope of recovering useful information about the [high-dimensional parameter] vector with limited samples” hastie2015statistical. Imposing some constraints, hopefully driven by economic reasoning (as we promote here), should help to extract relevant information in big data problems. As a consequence, the aOGL approach yields better performance in terms of covariate selection and estimated models without arbitrage (see our Monte Carlo results in Section (ref) and our empirical results in Section (ref)). On our data for US single stocks, more than half of the stocks require dynamics in their factor loadings, while penalization without (with) grouping yields to 100% (0%) of all estimated time-varying models violating the no-arbitrage restrictions. Besides, the aOGL approach yields better in-sample and out-of-sample predictive performance on an equally-weighted portfolio (see Sections (ref) and (ref)). On our data for US single stocks, prediction errors are located closer to zero and their scale is narrower.

LASSO type techniques have already been applied successfully to factor models in finance. bryzgalova2015spurious develops a shinkrage-based estimator that identifies the weak factors (i.e., factors that do not correlate with the assets) and ensures consistent and normality of the estimates of the risk premia. feng2020taming propose a model-selection method to evaluate the risk prices of observable factors. freyberger2020dissecting propose a nonparametric method to determine which firm characteristics provide incremental information for the cross section of expected excess returns. gu2020empirical use penalisation techniques for prediction purposes. Alternatively, Fan2022 develop a nonparametric methodology for estimating conditional asset pricing models using deep neural networks, by employing time-varying conditional information on alphas and betas carried by firm-specific characteristics. ACMV2022 propose a novel Bayesian approach to study time-series and cross-sectional effects in asset returns, when the true factor model and its underlying parameters are uncertain. They use macro predictors to model time-variation in the factor loadings and investigate potential mispricing. While their prior beliefs are weighted against mispricing, their analysis shows that time-varying mispricing appears with a large probability. CPZ2022 use deep neural networks to estimate a stochastic discount factor model for individual stock returns and exploit the fundamental no-arbitrage condition as criterion function, to construct the most informative test assets with an adversarial approach. Finally, let us mention that there is also work on inference for large dimensional models with observable and unobservable factors with high frequency data fan2016incorporating,ait2017using, pelger2019state, ait2020inference.

The outline of this paper is as follows. Section (ref) describes the conditional linear factor models with sparse time-varying coefficients, and how to implement the no-arbitrage restrictions in the specification of the random coefficient panel model. Section (ref) develops our penalized two-pass regression with time-varying factor loadings. The penalization in the first-pass (time-series) regressions of Section (ref) enforces sparsity for the time-variation drivers while also maintaining compatibility ex-ante with the no-arbitrage restrictions through building appropriate groups of coefficients. We explain in detail in Section (ref) why we prefer the aOGL method over the original Group-LASSO of yuan2006model for the first-pass regression. The second-pass (cross-sectional) regression of Section (ref) delivers risk premia estimates to predict equity excess returns. In Section (ref), we show asymptotic consistency of our penalised two-pass regression estimates under an adaptive estimation for the first-pass regression coefficients. Section (ref) reports our simulations results. Section (ref) gathers our empirical results. After describing our data on US single stocks in Section (ref), we present our empirical results on in-sample and out-of-sample prediction performances in Sections (ref) and (ref). We investigate 13 characteristics and 6 common instruments for the dynamics of factor loadings, and use the four-factor model of carhart1997persistence and the five-factor model of fama2015five. Section (ref) concludes. We list regularity conditions in Appendix (ref), the proofs of our theoretical results in Appendices (ref) and (ref), and a description on how to construct groups for the numerical optimisation in Appendix (ref).

Model specification

In this section, we consider a conditional linear factor model with time-varying coefficients as in \citetalias{gagliardini2016time} (see gagliardini2019estimation for a review). From their Assumptions APR.1, APR.2, and APR.3, the time-varying factor model for assets belonging to the continuum of assets $\gamma \in [0,1]$ is

where $R_{t} (\gamma)$ denotes the excess return on asset $\gamma$ at period $1, \ldots, T$, vector $f_{t} \in \mathbb{R}^K$ gathers the values of the factors at date $t$. From Assumption APR.1 of \citetalias{gagliardini2016time}, the intercept $a_{t}(\gamma) \in \real$ and factor loadings $b_{t}(\gamma) \in \real^K$ are $\mathcal{F}_{t-1}$-measurable, where the filtration process $\mathcal{F}_{t-1}$ is the information available to all investors at time $t-1$. The error terms have mean zero $\mathbb{E}[\varepsilon_{t}(\gamma) \vert \mathcal{F}_{t-1}] = 0$ and are uncorrelated with the factors conditionally on information $\mathcal{F}_{t-1}$, ${\operatorname*{Cov}(\varepsilon_{t}(\gamma), f_{t,k} \vert \mathcal{F}_{t-1})} = 0$, $k=1,...,K$. Assumption APR.2 of \citetalias{gagliardini2016time} gathers standard measurability conditions for a stochastic process, and requires that the process $\beta_t(\gamma) = (a_{t}(\gamma),b_{t}(\gamma)^\top)^{\top} \in \real^{K+1}$ is a bounded aggregate process as defined in AlNajjar_1995, as well as the nondegeneracy in the factor loadings across assets. Assumption APR.3 of \citetalias{gagliardini2016time} imposes an approximate factor structure in ((ref)) such that, for any sequence $\gamma_i \in [0,1], i = 1, \ldots, n$, with $\Sigma_{\varepsilon_{t},t,n} \in \real^{n\times n}$ being the conditional variance-covariance matrix of the vector $(\varepsilon_{t}(\gamma_1), \ldots ,\varepsilon_{t}(\gamma_n))^{\top}$ knowing $Z_{t-1}$, there exists a set such that $n^{-1}\operatorname{eig}_{\text{max}}(\Sigma_{\varepsilon_{t},t,n} ) \overset{\text{\tiny{$L^2$}}}{\longrightarrow}0$ as $n\to\infty$, where $\operatorname{eig}_{\text{max}}(\Sigma_{\varepsilon_{t},t,n})$ denotes the largest eigenvalue of $\Sigma_{\varepsilon_{t},t,n}$, and where $\overset{\text{\tiny{$L^2$}}}{\longrightarrow}$ denotes convergence in the $L^2$-norm. Under Assumptions APR.4 of \citetalias{gagliardini2016time}, the following asset pricing restriction holds:

for all $\gamma \in [0,1]$, at any date $t = 1, 2, \ldots$ where random vector $\nu_{t}\in\mathbb{R}^{K}$ is unique and is $\mathcal{F}_{t-1}$-measurable, which can also be written as

with $\lambda_{t}=\nu_{t}+\mathbb{E}[f_{t}\vert\mathcal{F}_{t-1}] \in \real^K$. Equation ((ref)) shows the link between expected excess returns and the product of the time-varying factor loadings and risk premia. Below, we rely on that link to predict excess returns. Assumption APR.4 of \citetalias{gagliardini2016time} excludes asymptotic arbitrage opportunity, such that there is no portfolio sequence with zero cost and positive payoff. The conditioning information $\mathcal{F}_{t-1}$ contains $Z_{\underline{t-1}}$ and $Z_{\underline{t-1}}(\gamma)$, where $Z_{t-1}\in\mathbb{R}^{p}$ is a vector of lagged instruments common to all stocks, $Z_{t-1}(\gamma)\in\mathbb{R}^{q}$, for $\gamma \in [0,1]$, is a vector of lagged characteristics specific to stock $\gamma$, and $Z_{\underline{t}}=\{Z_{t},Z_{t-1},...\}$ denotes the set of past realizations. Vector $Z_{t-1}$ may include past observations of the factors and some additional variables such as macroeconomic variables. Vector $Z_{t-1}(\gamma)$ may include past observations of firm characteristics and stock returns. We define the dynamics of the factor loadings $b_t(\gamma)$ as a sparse linear function of $Z_{t-1}$ Shanken_1990,Ferson_Harvey_1991 and $Z_{t-1}(\gamma)$ Avramov_Chordia_2006.

Moreover, we define the vector of risk premia as a sparse linear function of lagged instruments $Z_{t-1}$ Cochrane_1996,Jagannathan_Wang_1996 and specify the conditional expectation of the factor $\mathbb{E}\left[ f_t \vert \mathcal{F}_{t-1} \right]$ given the filtration process $\mathcal{F}_{t-1}$.

Assumptions (ref) and (ref) differ from Assumptions FS.1 and FS.2 of \citetalias{gagliardini2016time}. Indeed, we consider here the matrices $B(\gamma), C(\gamma), \Lambda_1$ and $F_1$ of coefficients as sparse, meaning that only a small fraction of the $Z_{t-1}$ or $Z_{t-1}(\gamma)$ for $\gamma \in [0,1]$ are useful to describe the dynamics of the factor loadings, risk premia, and conditional expectation of the factors. Building on the sampling scheme from Assumptions SC.1 and SC.2 of \citetalias{gagliardini2016time}, we define the indicator variable $I_t(\gamma)$, for all $\gamma \in [0,1]$, such that $I_t(\gamma) = 1$ if the return on asset $\gamma$ is observable at time $t$, and 0 if not. Assumption SC.1 ensures that $I_t(\gamma)$, $\varepsilon_t(\gamma)$ and variables in $\mathcal{F}_{t-1}$ are independent, while Assumption SC.2 ensures that the random variables $\gamma_{i}$, $i=1,...,n$, are i.i.d.\ indices, independent of $\varepsilon_{t}(\gamma)$, $I_{t}(\gamma)$, and $\mathcal{F}_{t-1}$. From the above sampling scheme, we can now use the following notation: $I_{i,t} = I_t(\gamma_i), R_{i,t} = R_t(\gamma_i),\beta_{i,t} = \beta_t(\gamma_i), \varepsilon_{i,t} = \varepsilon_{t}(\gamma_i), A_i = A(\gamma_i), B_i = B(\gamma_i), C_i = C(\gamma_i)$ and $Z_{i,t-1} = Z_{t-1}(\gamma_i) $ as well as $a_{i,t} = a_t(\gamma_i)$ and $b_{i,t} = b_t(\gamma_i)$. Hence, from Assumptions (ref) and (ref), we can express ((ref)) using the asset pricing restriction in ((ref)) as the following Data Generating Process (DGP):

We see that the first term $A_i^{\top} \left(\Lambda_0 - F_0\right)$ corresponds to the time-invariant part in the time-varying intercept $a_{i,t}$, while the term $A_i^{\top}f_{t}$ corresponds to the time-invariant part of the time-varying factor loadings $b_{i,t}$. To separate the time-invariant part from the time-varying part, we make the following assumption on the model specification.

Assumption (ref) ensures that the time-invariant part of a factor loading is always included in the model specification, so that we can distinguish a factor with a time-invariant loading from a factor with a time-varying loading for asset $i$. This assumption is key to analyze which instrument $Z_{t-1}$ and characteristic $Z_{i,t-1}$, if needed, drive the dynamics of the factor loadings $b_{i,t}$ for asset $i$, and impact on the prediction $\mathbb{E}[R_{i,t}\vert\mathcal{F}_{t-1}]$ via ((ref)). Since implementing a penalized two-pass regression given on ((ref)) is difficult (due to the quadratic form in lagged instruments $Z_{t-1}$ and $Z_{i,t-1}$), we redefine the regressors and coefficients, as a generic panel model. Beforehand, let us define the vector of lagged instruments including the intercept as $\tilde{Z}_{t-1} = (1, Z_{t-1}^{\top})^{\top} \in \real^{\tilde{p}}$, where $\tilde{p} = p+1$, and the new matrices $\breve{B_i} = [A_i | B_i ] \in \real^{K \times \tilde{p}} $ and ${\Lambda} - {F} = [(\Lambda_0 - F_0) | (\Lambda_1 - F_1) ] \in \real^{K \times \tilde{p}} $ that stack respectively column-wise the elements of $A_i$, $B_i$, and $(\Lambda_0 - F_0) ,(\Lambda_1 - F_1 )$. The linear transformed regressors are

where $d_2 = d_{21} + d_{22} = K\tilde{p} + Kq$, and

where $d_{1} = d_{11} + d_{12} =(\tilde{p}+1)\tilde{p}/2+\tilde{p}q$ and the symmetric matrix $X_{t}=(X_{t,k,l})_{k,l}\in\mathbb{R}^{\tilde{p}\times \tilde{p}}$ is such that $X_{t,k,l}=\tilde{Z}_{t-1,k}^{2}$, if $k=l$, and $X_{t,k,l}=2\tilde{Z}_{t-1,k}\tilde{Z}_{t-1,l}$, otherwise, for $k,l=1,\ldots,\tilde{p}$, where $\tilde{Z}_{t,k}$ denotes the $k$-th component of the vector $\tilde{Z}_{t}$. The vector-half operator $\operatorname*{vech}\left[\cdot\right]$ stacks the elements of the lower triangular part of a $\tilde{p}\times \tilde{p}$ matrix as a $\tilde{p}\left(\tilde{p}+1\right)/2 $ vector. The first element of $\operatorname*{vech}\left(X_{t}\right)$ is related to the time-invariant coefficients $A_i^{\top} \left(\Lambda_0 - F_0\right)$, whereas the elements $2, \ldots, \tilde{p}$ are related to $A_i^{\top} \left(\Lambda_1 - F_1\right) Z_{t-1} + Z_{t-1}^{\top} B_i^{\top} \left(\Lambda_0 - F_0\right)$. Through the above redefinitions of the regressor, we can write ((ref)) as

where $x_{i,t} = (x_{1,i,t}^{\top}, x_{2,i,t}^{\top})^{\top}$ is of dimension $d = d_1 + d_2$ and $\beta_i = (\beta_{1,i}^{\top}, \beta_{2,i}^{\top})^{\top}$ is defined as

and where $W_{\tilde{p},q}$ is the commutation matrix such that $\operatorname*{vec}[M^{\top}] = W_{\tilde{p},q} \operatorname*{vec}[M]$. Moreover, $D_{\tilde{p}}^{+}$ denotes the $((\tilde{p}+1)\tilde{p}/2+\tilde{p}q) \times \tilde{p}^2$ Moore-Penrose inverse of the duplication matrix $D_{\tilde{p}}$ such that $\operatorname*{vech}[M] = D_{\tilde{p}}^{+} \operatorname*{vec}[M]$, for any matrix $\tilde{p} \times \tilde{p}$ matrix $M$. The following section describes the selection and estimation part of the model.

Estimation and selection

This section implements the two-pass regression of jensen1972capital and fama1973risk, while selecting the contributing variables in the time-varying factor loadings. The penalized first-pass (time-series) regression selects and estimates the non-zero coefficients $\beta_i$ for $i = 1, \ldots, n$, ensuring a model specification compatible ex-ante with the no-arbitrage restrictions through the aOGL approach of percival2012theoretical. The second-pass regression relies on the Weighted Least-Square (WLS) estimator of \citetalias{gagliardini2016time} to estimate the vector $\nu$, and takes the adaptive LASSO (aLASSO) estimator of zou2006adaptive to select and estimate the matrix $F$ of coefficients of the conditional expectation of the factors.

First-pass regression

The goal of the penalized first-pass regression is to select and estimate the factor loadings for each asset $i = 1, \ldots, n$, while keeping their respective time-invariant contribution fully specified as described in Assumption (ref). Moreover, it aims at selecting variables ensuring a proper model specification consistent ex-ante with the no-arbitrage restrictions for each stock. A possible solution to ensure that these restrictions are satisfied while allowing to select variables in the first-pass regression is to consider a LASSO-type estimator based on appropriate predefined sets of indices corresponding to groups of variables. We define $\groupset \subset \mathcal{P}(\{1, \ldots, d\})$ as the set of indices corresponding to all possible (potentially overlapping) groups in line with the no-arbitrage restrictions, where $\mathcal{P}(\{1, \ldots, d\})$ denotes the power set of $\{1, \ldots, d\}$. Moreover, we let $g \in \groupset$ denote a possible group and we require that the indices associated to all covariates belong to at least one group. Under the framework discussed in the previous sections, we define below the restrictions on $\groupset$ such that a model selection procedure based on $\groupset$ satisfies ex-ante the no-arbitrage restrictions by construction. \setcounter{Restriction}{0}

These restrictions ensure that Assumption (ref) is satisfied and that a model selection procedure guarantees that the instrument $\tilde{Z}_{t-1,l}$ or characteristic $Z_{i,t-1,m}$ exist in either both $x_{1,i,t}$ and $x_{2,i,t}$, or neither. More specifically, Restriction (ref) is related to Assumption (ref), which requires the coefficients in $\beta_i$ related to the time-invariant contribution to be always included in the selected model. Restriction (ref) is related to Assumption (ref) and Assumption (ref). Under the DGP in ((ref)), and from the definition of $\operatorname*{vech}(X_t)$, we can see that the off-diagonal of $X_t$ in $\operatorname*{vech}[X_t]$ cannot be assigned to any groups. We cannot assign $2\tilde{Z}_{t-1,s}\tilde{Z}_{t-1,l}$ to a group a priori, since its contribution can come from either the specification in Assumption (ref) or (ref). Restriction (ref) reflects this point, and imposes no specific group-structure to those covariates which are penalized individually. Restrictions (ref) and (ref) are critical in the model building. They constrain the set of possible models only to those compatible with the no-arbitrage restrictions, so that we do not introduce arbitrage ex-ante in the model specified in ((ref)). We want to avoid that the no-arbitrage restriction $a_{i,t} = b_{i,t}^{\top} \nu_t$ is violated by construction ex-ante in the specification.

To satisfy the above restrictions, the Group-LASSO of yuan2006model constrains the set of possible models. For its implementation, we need to create a group with all scaled factors and their corresponding terms in the intercept, hence it implies that we select either all scaled factors (keep the group) or none of them (delete the group). To illlustrate this point, let us consider the following simple case with one common instrument, say inflation, and the Fama-French five-factor model fama2015five. The Group-LASSO would force us to select either all scaled factors (product between lagged inflation and the factors), or none of them. It removes the possibility that only a subset of them is relevant; for example, only the product of inflation and the market factor matters for the dynamics of excess returns. Besides, we could think of using multiple groups, each one containing one scaled factor and its associated instrument. jacob2009group investigate such a proposal and show that this approach is not appropriate as the Group-LASSO removes all groups if at least one of those groups is not selected.

To tackle this problem, jacob2009group propose the OGL, or latent Group-LASSO. They introduce the latent variables ${v}_{\group} \in \mathcal{V}_g = \{ x \in \real^d | \operatorname*{supp}(x) = g \}$, for $g \in \groupset$ and where $\operatorname*{supp}(x)$ denotes the support of $x$, i.e., the set of indices $i \in \{1, \ldots, d\}$ such that $x_i \neq 0$. Moreover, we define ${v_g} = (v_{g_1}^\top, \ldots, v_{g_{J}}^\top)^\top \in \real^d, \mathcal{V}(\beta) = \{v_g : g\in \groupset\}$, s.t.\ $\beta = \sum_{g\in \groupset} {v}_{g}$, and $J = |\groupset|$, $|\cdot|$ denotes the cardinality of a set and $g_j, \, j = 1,\ldots, J$, denotes the $j$-th element of $\groupset$. Hence, the OGL estimator is the solution of the following optimization problem:

with the penalty term $\|\beta_i\|_{2,1,\groupset}$ defined as

where $\|\cdot\|$ denotes the $l_2$-norm. In this work, we consider the adaptive version of OGL (aOGL) studied by percival2012theoretical, for which the estimator is described as follow:

where $\delta_{\group} \geq 0$ denotes the data-dependent (adaptive) weight associated to group $g$, and $\delta \geq 0$ corresponds to the overall amount of shrinkage. There are different strategies available in the literature for the Group LASSO and OGL to get estimator consistency and support selection consistency. They are based on the irrepresentable condition bach2008consistency,jacob2009group, adaptive shrinkage nardi2008asymptotic,percival2012theoretical and group sparsity lounici2011oracle. We choose adaptive shrinkage since it simplifies the presentation and derivation of our asymptotic results in a random design setting. Since our goal is to shrink toward the model that includes only the time-invariant contribution of the covariates, the weight associated with the first element of $\delta_{\group}$ is equal to zero. The penalty term in ((ref)) leads to a solution which is a union of the groups due to the latent variables ${v}_{g}$. One strategy to solve the minimization problem given in ((ref)) and ((ref)) is the duplication of covariates put forward in jacob2009group, that we adapt to our setting. In line with Restrictions (ref) to (ref), we consider 4 different group types. The first group includes the time-invariant intercept and time-invariant factors, and is not penalised. The second set of groups contains the covariates related to Restriction (ref), which are penalized individually. The next two sets of groups consider Restrictions (ref) and (ref). They respectively group the terms in $Z^2_{t-1}$ and $Z_{i,t-1}$ from $x_{1,i,t}$ with their corresponding scaled factors in $x_{2,i,t}$. The columns of the initial vector with the elements indexed by the group $g$, which need to be duplicated, create a new vector of duplicated regressors. Then, we can solve the optimization problem in ((ref)) considering the duplicated regressors (instead of the initial ones), using the existing standard algorithm for the Group-LASSO. Appendix (ref) describes in detail how to construct those groups complying with the no-arbitrage restrictions ex-ante, and yielding the full vector of duplicated regressors used in the numerical optimisation.

Let us now compare the number of possible models under aOGL and aLASSO methods. For the aOGL approach, we can associate a model to every subset of $\mathcal{G}$. Indeed, consider $\mathcal{W} \subseteq \groupset$, then this subset is associated to the set $S_{\mathcal{W}} = \bigcup_{l = 1}^{|\mathcal{W}|} \mathcal{W}_l$ of indices. It allows us to enumerate the number $2^{J-1}$ of possible models under appropriate grouping. That number is typically much lower in empirical applications than the number $2^{d-n_1}$ of possible models with a LASSO penalization, where $n_1$ is the number of covariates associated to the time-invariant contribution group. We get the ratio $2^{J-1}/2^{d-n_1} = 2^{-(pq+p+q)}$, and we can see that, for large $p$ and $q$, the aLASSO method examines many more possibilities. Besides, from Assumption (ref), we have $\min(p,q) \geq 1$, and deduce the upper bound:

To further illustrate the grouping structure and the importance of Restrictions (ref) to (ref), let us consider the following simple two-factor model with a single common instrument and a single characteristic. Here, we have $K=2$, $\tilde{p} = 2$, and $q = 1$, with $\tilde{Z}_{t-1} = (1,Z_{t-1})^{\top} \in \real^2$, so that the regressors $x_{i,t} = (x_{1,i,t}^{\top},x_{2,i,t}^{\top})^{\top}$ become

and

with their respective coefficients $\beta_{1,i} = (\beta_{1,i,1},\beta_{1,i,2},\beta_{1,i,3},\beta_{1,i,4},\beta_{1,i,5})^{\top}$ and $\beta_{2,i} = (\beta_{2,i,1},\beta_{2,i,2},\beta_{2,i,3},\beta_{2,i,4},\beta_{2,i,5},\beta_{2,i,6})^{\top}$.

\sloppy From the definition of grouping structure in Apprendix (ref), we construct the set of six groups made of the covariates: \break $(x_{1,i,t,1}, x_{2,i,t,1}, x_{2,i,t,3})^{\top}$ for the time-invariant contribution, $(x_{1,i,t,2})$ for the covariate associtated to Restriction (ref), $(x_{1,i,t,3}, x_{2,i,t,2})^{\top}$ and $(x_{1,i,t,3}, x_{2,i,t,4})^{\top}$ grouping the covariates in $\tilde{Z}_{t-1}$, and finally $(x_{1,i,t,4}, x_{1,i,t,5}, x_{2,i,t,5})^{\top}$ and $ (x_{1,i,t,4}, x_{1,i,t,5}, x_{2,i,t,6})^{\top}$ grouping the covariates in $\tilde{Z}_{i,t-1}$. Stacking those vectors row-wise in a single column defines the full vector of duplicated covariates for the numerical optimisation in the aOGL estimation. Besides, we can use this simple example to illustrate two possible manners to introduce ex-ante arbitrage through careless modeling. Removing the covariates $x_{2,i,t,2} = Z_{t-1} f_{t,1}$ and $x_{2,i,t,4} = Z_{t-1} f_{t,2}$ from the full model might introduce ex-ante arbitrage through $x_{1,i,t,3} = Z_{t-1}^2$ since we miss its associated scaled factors in $x_{2,i,t}$. Here, the coefficient associated with $x_{1,i,t,3}$ might be shrunk to zero by the aLASSO estimator, avoiding ex-post a model with arbitrage. On the contrary, removing the quadratic term $x_{1,i,t,3}$, while keeping its corresponding scaled factors $x_{2,i,t,2}$ and $x_{2,i,t,4}$, introduces ex-ante arbitrage in the model by construction, since we cannot estimate the coefficient of $x_{1,i,t,3}$, when that covariate is absent from the model.

Table (ref) lists the set $\model = \{\model_1, \ldots, \model_{32}\}$ of possible models that respect Restrictions (ref) to (ref) with $\model_1$ being the model with the time-invariant contribution only (Assumption (ref)). The aOGL method gives $2^5$ possible models. It is considerably smaller than the $2^{8} = 256$ possible models under the aLASSO method. Here, we reach the upper bound ((ref)) since $p = q = 1$. We can see that our regularization approach restricts the space of searched models, even in this simple time-varying setting, and hence permits a sound exploration of the possible models consistent with finance theory. Moreover, the two specifications with arbitrage described in the above lines are not in the set $\model$ of models induced by the grouping structure of the aOGL approach, strengthening conducive arguments for our proposed method.

Having showed the advantages of the aOGL in terms of model building, we now state the asymptotic result of the first-pass regression. Beforehand, we introduce some notations from percival2012theoretical. Let us define the two sets of indices \break $H_i=\left\{l \in \{1, \ldots, d\}: \beta_{i,l} \neq 0\right\}$, $H^c_i= \{l \in \{1, \ldots, d\}: \beta_{i,l} = 0\}$, corresponding to the sets of non-zero and zero true coefficient $\beta_i$. Moreover, we take

the sets of groups in which the indices are respectively all non-zero, all zero and a mix of zero and non-zero in $\beta_i$. To investigate the asymptotic properties of the estimator in ((ref)), we make the following assumptions:

Assumption (ref) is a usual assumption for the standard OLS solution to be consistent (see Assumption B.1 in \citetalias{gagliardini2016time}), while Assumption (ref) allows for or a martingale difference sequence for the error terms (see Assumption A.1 in \citetalias{gagliardini2016time}). Assumption (ref) is discussed in percival2012theoretical and addresses the uniqueness of decomposition of $\beta_i$. We now state the main result for the first-pass regression, which corresponds to Theorem 2 derived by percival2012theoretical in the fixed design framework.

For the above result to hold, the vector $\beta_i^{\text{init}}$ needs to be $\sqrt{T_i}$-consistent. More specifically, $\beta_i^{\text{init}}$ is any $a_{T_i}$-consistent estimator where $a_{T_i} \to \infty$, and $a_{T_i}^{\check{\gamma}} \sqrt{T_i} \delta \to \infty$. Moreover, in the context of the aOGL, the decomposition $\{v^{\text{init}}_{i,g}\}$ must be unique. Lemma 4 in percival2012theoretical shows $\sqrt{T_i}$-consistency of the $\{v^{\text{OLS}}_{i,g}\}$, which is a example of a potential solution for $\{v^{\text{init}}_{i,g}\}$ in the case of fixed covariates. In our framework, with the uniqueness assumption of the decomposition, we can use the ridge regression estimator as $\{v^{\text{init}}_{i,g}\}$. The distributional result of Lemma (ref) is key in deriving the asymptotic properties of the second-pass regression discussed in the next section.

To control for short sample size, and potentially numerical instability on the inversion of matrix $\hat{Q}_{x,i}$, we consider the trimming device defined in \citetalias{gagliardini2016time}, such that $\boldsymbol{1}_{i}^{\chi}=\boldsymbol{1}\{ CN(\hat{Q}_{x,i})\leq\chi_{1,T},\tau_{i,T}\leq\chi_{2,T}\} $, where $ CN(\hat{Q}_{x,i})=\sqrt{\operatorname{eig}_{\max}(\hat{Q}_{x,i})/\operatorname{eig}_{\min}(\hat{Q}_{x,i})}$ is the condition number of the matrix $\hat{Q}_{x,i}$, $\operatorname{eig}_{\min}(\cdot)$ denotes the minimum eigenvalue, and $\tau_{i,T}=T/T_{i}$. The first trimming based on $CN(\hat{Q}_{x,i})\leq\chi_{1,T}$ selects the assets for which the time-series regression is not badly conditioned, while the second trimming based on $\tau_{i,T}\leq\chi_{2,T}$ keeps only the assets for which samples are not too short.

Second-pass regression

The second-pass regression aims at computing the cross-sectional estimator of $\nu$. For that purpose, we implement the WLS estimator of \citetalias{gagliardini2016time}, while accounting for the sparse model specification in the first-pass regression for all $i = 1, \ldots, n$. For that purpose, we introduce the indicator vector $ \mathbf{1}_{\beta_i} \in \natural^d$, such that $\mathbf{1}_{\beta_{i,l}} = 1$ if $\beta_{i,l} \neq 0$, and $0$ otherwise, for $l = 1, \ldots, d$, that we decompose in the following manner: $ \mathbf{1}_{\beta_i} = (\mathbf{1}_{{\beta}_{11,i}}^{\top}, \mathbf{1}_{{\beta}_{12,i}}^{\top}, \mathbf{1}_{{\beta}_{21,i}}^{\top}, \mathbf{1}_{{\beta}_{22,i}}^{\top})^{\top}, $ where $\mathbf{1}_{{\beta}_{11,i}} \in \natural^{d_{11}}$, $\mathbf{1}_{{\beta}_{12,i}} \in \natural^{d_{12}}$, $\mathbf{1}_{{\beta}_{21,i}} \in \natural^{d_{21}}$ and $\mathbf{1}_{{\beta}_{22,i}} \in \natural^{d_{22}}$. To implement the WLS estimator for the vector $\nu$, we need to account for the different number of regressors selected through the aOGL approach. Hence, in the same spirit as in chaieb2020factors, we introduce the following selection matrices that help us transforming the $x_{i,t}$ into their sparse counterparts. The matrices $\tilde{D}_i$ and $\tilde{E}_i$ are the $d_{11} \times d_{11,i}$ and $d_{12} \times d_{12,i}$ such that columns with all zeros have been removed in $\operatorname*{diag}[\mathbf{1}_{{\beta}_{11,i}}]$ and $\operatorname*{diag}[\mathbf{1}_{{\beta}_{12,i}}]$. Similarly, the matrices $\tilde{B}_i$ and $\tilde{C}_i$ are the $d_{21,i} \times d_{21}$ and $d_{22,i} \times d_{22}$ matrices such that rows with all zeros have been removed in $\operatorname*{diag}[\mathbf{1}_{{\beta}_{21,i}}]$ and $\operatorname*{diag}[\mathbf{1}_{{\beta}_{22,i}}]$. Moreover, we define $x_{H_i,i,t}$ as the vector of regressors indexed by $H_i$ after the selection of the first pass.

Based on the selection matrices $\tilde{D}_i,\tilde{E}_i,\tilde{B}_i$, and $\tilde{C}_i$ , we rewrite the parameter restriction in ((ref)) such that

where $N_{\tilde{p}}$ is defined in ((ref)), yielding the asset pricing restrictions expressed in the newly defined ${\beta}_{1,i}$ and ${\beta}_{3,i}$ as $ {\beta}_{1,i} = {\beta}_{3,i} \nu,$ $ \nu=\operatorname*{vec}[\Lambda^{\top}-F^{\top}]. $ We obtain ${\beta}_{3,i}$ from the following identity,

We can now implement the following second-pass regression WLS estimator

where $\hat{\nu}$ denotes the estimator of $\nu$, $\hat{Q}_{\beta_{3}}=\frac{1}{n}\sum_{i}\hat{\beta}_{3,i}^{\top}\hat{w}_{i}\hat{\beta}_{3,i}$ , and weights are estimates of $w_{i}=\boldsymbol{1}_{i}^{\chi}\left(\operatorname*{diag}\left[v_{i}\right]\right)^{-1}$. Moreover, the $v_{i}$ are the asymptotic variances of the standardized errors $\sqrt{T}(\hat{\beta}_{1,i}-\hat{\beta}_{3,i}\nu)$ in the cross-sectional regression for large $T$ such that $ v_{i}=\tau_{i}C_{\nu,1,i}^{\top}Q_{H_i,{x},i}^{-1}S_{ii}Q_{H_i,{x},i}^{-1}C_{\nu,1,i}$, where $ S_{ii}=\operatorname*{plim}_{T\to\infty}\frac{1}{T}\sum_{t}\sigma_{i}^2{x}_{H_i,i,t}{x}_{H_i,i,t}^{\top}$ and $C_{\nu,1,i}$ $=$ $(E_{1,i}^{\top}-(I_{d_{1,i}}\otimes\nu^{\top})J_{a,i}E_{2,i}^{\top})^{\top}$, $E_{1,i}=(I_{d_{1,i}},0_{d_{1,i}\times d_{2,i}})^{\top}$, $E_{2,i}=(0_{d_{2,i}\times d_{1,i}},I_{d_{2,i}})^{\top}$. We use the estimates $ \hat{v}_{i}=\tau_{i,T}C_{\hat{\nu}_{1}}^{\top}\hat{Q}_{H_i,{x},i}^{-1}\hat{S}_{ii}\hat{Q}_{H_i,{x},i}^{-1}C_{\hat{\nu}_{1}}$, where $ \hat{S}_{ii}=\frac{1}{T_{i}}\sum_{t}I_{i,t}\hat{\varepsilon}_{i,t}^{2}{x}_{H_i,i,t}{x}_{H_i,i,t}^{\top}$, $\hat{\varepsilon}_{i,t}=R_{i,t}-\hat{\beta}_{i}^{\top}{x}_{H_i,i,t}$ together with $ C_{\hat{\nu},1,i}=(E_{1,i}^{\top}-\left(I_{d_{1,i}}\otimes\hat{\nu}_{1,i}^{\top}\right)J_{a,i}E_{2,i}^{\top})^{\top}$. To estimate $C_{\nu,1,i}$, we use the OLS estimator given by $ \hat{\nu}_{1,i}=(\sum_{i}\boldsymbol{1}_{i}^{\chi}\hat{\beta}_{3,i}^{\top}\hat{\beta}_{3,i})^{-1}\sum_{i}\boldsymbol{1}_{i}^{\chi}\hat{\beta}_{3,i}^{\top}\hat{\beta}_{1,i}$. We estimates the weights with $\hat{w}_{i}=\boldsymbol{1}_{i}^{\chi}\left(\operatorname*{diag}\left[\hat{v}_{i}\right]\right)^{-1}$.

To study the asymptotic properties of the estimator $\hat{\nu}$, we consider the following assumption on the size of the cross-section $n$.

Assumption (ref) puts a bound on the growth of the cross-section such that it does not grow faster that some power of the sample size $T$. In Proposition (ref), we provide the consistency result for the estimator $\hat{\nu}$.

Assumptions (ref) to (ref) are discussed in Appendix (ref). This asymptotic property of $\hat{\nu}$ is studied under the double asymptotics $n,T \to \infty$ in \citetalias{gagliardini2016time}. They show consistency of $\hat{\nu}$ under a full representation of $\beta_i$, while we assume a sparse representation of $\beta_i$. Hence, our result differs in that respect.

Let us now recover the sparse structure of the conditional expectation of the factors under Assumption (ref). For that purpose, we consider the aLASSO estimator of zou2006adaptive to select and estimate the matrix $F$ of coefficients. We solve the following minimization problem for all factor $f_{k,t},k = 1, \ldots, K$, such that the estimator of the $k$-th row of the matrix $F$ is given by:

where $\delta$ accounts for the overall amount of shrinkage as in ((ref)), and $\hat{w_j}$ are data dependent weights. Typically, the weights are defined as $\hat w_j = \nicefrac{1}{|\hat{F}_{1,k,j}^{\text{OLS}}|^{\check{\gamma}}}$ for $\check \gamma >0$, where $\hat{F}_{1,k,j}^{\text{OLS}}$ are the OLS estimates of $F_{1,k,j}$, the true values in the vector parameter $F_{1,k}$. The estimate $\hat{F}$ stacks row-wise the elements of $(\hat{F}_{0,k}, \hat{F}_{1,k})$ obtained from ((ref)). Under Assumption (ref), no amount of shrinkage is applied to $F_0$ in $F$, to always keep the time-invariant contribution in the model. We get the final estimates of the sparse matrix $\Lambda$ from the relationship $ \operatorname*{vec}[\hat{\Lambda}^{\top}]=\hat{\nu}+\operatorname*{vec}[\hat{F}^{\top}]$, which yields $\hat{\lambda}_{t}=\hat{\Lambda}Z_{t-1}$. To derive the asymptotic consistency of $\hat{\Lambda}$, we rely on Proposition (ref) for the estimator $\hat{\nu}$. Let us consider the following assumption:

Assumptions (ref) is a standard regularity assumption on the design matrix for linear regression model, in order to obtain a unique solution for $(F_{0,k},F_{1,k})$. Under the above Assumption (ref), and Proposition (ref), the following proposition gives the consistency result for the estimator $\hat{\Lambda}$.

Proof of Proposition (ref) is direct since from the definition of $\hat{\Lambda}$, $\Vert \operatorname*{vec} [\hat{\Lambda}^{\top} - \Lambda^{\top}]\Vert \leq \Vert \hat{\nu} - \nu \Vert + \Vert \operatorname*{vec}[\hat{F}^{\top} - F^{\top}]\Vert$. From Proposition (ref), we know that $\Vert \hat{\nu} - \nu \Vert = o_p(1)$. Moreover, the aLASSO estimator in ((ref)) is a special case of the estimator in ((ref)), where each group is a singleton. Hence, considering Assumptions (ref) and (ref) which are the counterpart of Assumptions (ref) and (ref) respectively, we get that the result of Lemma (ref) applies to ((ref)). Hence, $\Vert \operatorname*{vec}[\hat{F}^{\top} - F^{\top}]\Vert = o_p(1)$. Therefore, we get consistency of $\hat{\lambda}_t$, $\sup_t\Vert \hat{\lambda}_t - \lambda_t \Vert = o_p(1)$, under Assumptions (ref) and (ref).

Simulation study

In this section, we study how the selection and estimation procedures of Section (ref) perform in finite samples. This first simulation study aims at investigating the prediction and selection performance of the aOGL method and at comparing it with the aLASSO method in a very sparse environment (Assumptions (ref) and (ref)). To that purpose, we simulate $500$ replicates from the DGP in ((ref)) for a (randomly drawn) single asset $i$ with sample size $T_i = 500$. We split that full sample in a training subsample and a testing subsample of $450$ and $50$ observations. The testing set is used for out-of-sample prediction performance assessment, where we compare the realized excess returns $R_{i,t}$ with their predictions $\hat{R}_{i,t} = \hat b_{i,t}^{\top} \hat \lambda_t$ under the model estimated on the training set. Errors in ((ref)) are i.i.d.\ such that $\varepsilon_{i,t}\sim\mathcal{N}(0,\sigma^2)$, where $\sigma = 0.09$. We match the model specification described in our empirical study (Section (ref)) for the common instruments $Z_{t-1} \in \real^6$ and stock-specific instruments $Z_{i,t-1}\in \real^{13}$. For the factors, we use the Fama-French five-factor model fama2015five described in the next section, namely we condition w.r.t.\ the values $f_t$ observed in our empirical study for the five factors. We also condition w.r.t.\ the observed $Z_{t-1}$ and $Z_{i,t}$ for asset $i$ of our empirical study. We only draw the error terms as in a parametric bootstrap.

In accordance with sparsity in Assumptions (ref) and (ref) and non-sparse time-invariant contribution in Assumption (ref), we set the matrices $A_i$, $B_i$, and $C_i$ according to their values for asset $i$ in the empirical study, with one non-zero element in $B_i$ and two non-zero element for $C_i$. We keep the vector $A_i$ full. We set the corresponding $a_{i,t}$ in order to avoid ex-ante arbitrage. Since we take very sparse matrices $B_i$ and $C_i$, we can view the simulation study as conservative for selection performance assessment (type of worst-case scenario). The resulting $\beta_i$ has 28 non-zero coefficients (including the 6 coefficients induced by the non-sparse time-invariant contribution) over a total of 219 coefficients. The matrices $F$ and $\Lambda$ are simply set to zero since they do not concern the aOGL estimator.

The selection and prediction performance is measured through the average Root Mean Squared Prediction Error (Av($ \text{RMSPE}_{R} $)), the average Root Mean Squared Error for parameter $\beta_i$ (Av($\text{RMSE}_{{\beta}}$)), the proportion of times the model introduces arbitrage (Arb.\ ($\%$)), the average number of selected true non-zero coefficients (True+), and average number of regressors in the selected model (NbReg). Table (ref) summarizes the results. The aOGL method makes a better job at predicting out-of-sample with a reduction of 1.7% w.r.t.\ the aLASSO method. The improvement in the average of RMSE for $\beta_i$ is 109%. The standard errors are also much lower (reduction of 9.1% and 91.0% for the Av($ \text{RMSPE}_{R}$) and Av($\text{RMSE}_{{\beta}}$)). Contrary to the aLASSO method, for which $98.2\%$ of estimated models exhibit arbitrage, the aOGL method selects only models without introducing ex-ante arbitrage by construction. Since we face less than $100\%$ for the aLASSO method, it sometimes shrinks adequately to zero the coefficients that should be. The aOGL method is able to recover in average the 11 true non-zero coefficients (11.26) while the aLASSO method struggles (7.37). The aOGL method is also more parsimonious than the aLASSO method in terms of selected regressors (average of 14.75 versus 16.05).

Our second simulation set-up focuses on the out-of-sample prediction performance of the aOGL method in a setting close to our empirical study of Section (ref). We use a training sample to estimate the model and a testing sample to gauge its out-of-sample prediction performance on an equally-weighted portfolio. We consider the same model specification in terms of $f_t$, $Z_{t-1}$ and $Z_{i,t-1}$ as in the first study and implement the following procedure. We sample randomly a subset of $n = 500$ assets from Section (ref) (training sample), while keeping the same proportion of time-invariant models as in Table (ref). From each asset $i$ in this subset, we simulate $T_i$ observations from $R_{i,t} = a_{i,t} + b_{i,t}^{\top} f_t + \varepsilon_{i,t}$ with the coefficients $a_{i,t}$ and $b_{i,t}$ chosen as their aOGL corresponding values for stock $i$. The $500 \times 1$ error vector $\varepsilon_{t}$ at date $t$ is Gaussian with mean zero and block-diagonal correlation matrix with 10 blocks of equal size $50$, where, within each block matrix, the correlation between $\varepsilon_{k,t}$ and $\varepsilon_{l,t}$ is set to $\text{corr}(\varepsilon_{k,t},\varepsilon_{l,t}) = 0.25^{|k-l|}$, $k,l= 1,..., 50, l \neq k$. The variance of each error $\varepsilon_{i,t}$ is set equal to 0.05. From those 500 simulated paths, we implement the aOGL estimation procedure of Section (ref), and compare it with the same procedure, but using the aLASSO estimator instead of the aOGL estimator to select the covariates in ((ref)).

To evaluate the out-of-sample prediction performance, we simulate one new cross-sectional sample (testing sample) from the time-varying factor model for the 500 assets and each date $t$ and compute the prediction $\hat{R}_{i,t} = \hat b_{i,t}^{\top} \hat \lambda_t$ for the 500 stocks and each date $t$ based on the estimator computed before through the aOGL and aLASSO methods. We finally compute the out-of-sample Prediction Error (PE) for an equally-weighted portfolio through the difference between the new simulated $\frac{1}{500} \sum_i R_{i,t}$ and its predicted value $\frac{1}{500} \sum_i \hat R_{i,t}$. We compute the Root Mean Squared Prediction Error (RMSPE), and the Mean Absolute Prediction Error (MAPE) over the vector gathering the PE at each out-of-sample date. We repeat this procedure 100 times to get an average and to compute a standard error. They are reported in Table (ref). We can see that the aOGL method is much better at out-of-sample predicting excess returns of an equally-weighted portfolio both in terms of average of MAPE (reduction by 14%) but also in terms of variability as measured by the standard errors (reduction by 84%). The empirical distribution of the prediction errors is given in Figure (ref). We can see that the aOGL method is centered closer to zero and with a lower dispersion when compared to the aLASSO method. Those second simulation results again point in favor of our advocated estimation method.

Empirical results

This section investigates the predictive capacity of the aOGL estimator and compares it with the aLASSO estimator. We also consider a pure time-invariant model, and a (hybrid) model with constant $\nu$ and time-varying risk premia. We use the aLASSO estimator to gauge the added value of incorporating the no-arbitrage restrictions in the penalisation approach and the time-invariant models to gauge the added value of allowing for full time-variation. The latter comparison checks that, when it comes to return prediction, the complicated model does not necessarily outperform because of potential overfitting.

Data description

We extract the stock returns from the CRSP database for US common stocks listed on the NYSE, AMEX, and NASDAQ, and remove stocks with prices below 5 USD. We exclude financial firms (Standard Industrial Classification Codes between 6000 and 6999). The firm characteristics come from COMPUSTAT. The sample begins in July 1963 and ends in December 2019. It gives us $T = 678$ monthly observations. We proxy the risk-free rate with the 1-month T-bill rate.

From freyberger2020dissecting, we consider the following $q=13$ firm level characteristics $Z_{i,t-1}$: change in share outstanding ($\Delta$ shrout), log change in the split adjusted shares outstanding ($\Delta$ so), growth rate in total assets (Inv), size (LME), last month volume over shares outstanding (lturnover), adjusted profit margin (PM), momentum and intermediate momentum ($r_{12,2}$ and $r_{12,7}$), short-term reversal ($r_{2,1}$), closeness to 52-week high (Rel_to_high), the ratio of market value of equity plus long-term debt minus total assets to Cash and Short-Term Investments (ROC), standard unexplained volume (SUV), and total volume (Tot_vol). We refer to freyberger2020dissecting for a detailed description of those characteristics. We only retain stocks for which all 13 characteristics are non-missing. It produces a sample of $n = 6874$. For each $Z_{i,t-1}$, we follow freyberger2020dissecting and compute the cross-sectional rank at each time $t-1$ for all observations chaieb2020factors. For the common instruments $Z_{t-1}$, we consider the $p=6$ following variables: dividend yield (dp), net equity expansion (ntis), inflation (infl), stock variance (svar), default spread (def_spread), and the term-spread (term_spread). For each $Z_{t-1}$, we center and standardize all observations.

We consider the two following sets of factors $f_t$. The first set is the four-factor model of carhart1997persistence, such that $f_t = (f_{m,t}, f_{hml,t}, f_{smb,t}, f_{mom,t})^{\top}$, where $f_{m,t}$ is the month $t$ market excess return over the risk free rate, $f_{hml,t}$, $f_{smb,t}$, $f_{mom,t}$ are respectively the month $t$ returns on zero investment factor-mimicking portfolio for size, book-to-market, and momentum. Our second set of factors considers the profitability factor $f_{rmw,t}$ and the investment factor $f_{cma,t}$ as in the five-factor model of fama2015five, such that $f_t = (f_{m,t}, f_{hml,t}, f_{smb,t}, f_{rmw,t},f_{cma,t})^{\top}$. Our choice for a parsimonious specification in the factor space is justified by our goal of studying the selection of common and idiosyncratic instruments $Z_{t-1}$ and $Z_{i,t-1}$ that have impacts on the dynamics of the $a_{i,t}$, $b_{i,t}$, and $\lambda_t$. gagliardini2019diagnostic and gagliardini2019estimation also report evidence that those factors with time-varying loadings are rich enough to achieve a weak cross-sectional dependence in the error terms, namely there are no remaining omitted factors in the error terms.

In-sample prediction performance and selection results

In this section, we investigate the selection results from the first-pass penalized regression. We compare the fit of the penalized two-pass procedure with aOGL described in Section (ref) to the aLASSO estimator, where we select the $x_{i,t}$ and estimate their coefficients in the first-pass regression with the aLASSO estimator of zou2006adaptive and fit the WLS estimator for the $\nu$ described in Section (ref). We compute the estimator $\hat{F}$ as in ((ref)). The horse race starts from the same set of initial data described in the previous section, and the comparison is thus made on the same initial full information. From the characteristics and common instruments outlined in Section (ref), under the Carhart four-factor model, we have $d=5$ for the time-invariant model and $d=199$ for the time-varying model. Regarding the five-factor model of fama2015five, we have $d=6$ and $d=219$ for the unconditional and conditional specifications. The number of possible models under the aLASSO method is $2^{194}$ ($2^{213}$) with $K=4$ ($K=5$), while the number of possible models under the aOGL method is $2^{97}$ ($2^{116}$), which gives the ratio $2^{-97}$, a much lower value than the upper bound $1/8$ in ((ref)).

We choose the regularisation parameter in a data dependent way to minimize the Akaike Information Criterion (AIC) for both aOGL and aLASSO estimator. As advocated in greene2008econometric, we use $\chi_{1,T} = 15$, and require at least 5 years of data such that $\chi_{2,T} = 678/60$. Because of the trimming, we do not keep the same set of stocks for each method and each model. Indeed, due to the different models induced by the first pass for each stock $i$, the trimming device $\boldsymbol{1}\{CN(\hat{Q}_{\check{x},i})\leq\chi_{1,T}\},$ yields a different set of stocks for each method. Since we do not wish to introduce multicolinearity in the second-pass regression, we choose to stick with different sets for each method. For the aOGL estimator, the aLASSO estimator, and the time-invariant estimator, we end up with 4412, 2225, 4879 for the four-factor model, and 4441, 2097, 4879 for the five-factor model. We can observe that the trimming device for the aLASSO method is more binding. As seen in the simulation results in Section (ref) and in Table (ref), the aLASSO method tends to include more variables, and, as a consequence, increase its associated condition number. Table (ref) reports the percentage (TI ($\%$)) of estimated models shrunk towards the time-invariant models. For those estimates, we only select the single group corresponding to Restriction (ref) related to Assumption (ref). Around two thirds of the stocks require dynamics in their factor loadings. This new empirical result based on a penalization approach illustrates the relevance of allowing for potential time-variation in modelling excess returns of individual stocks with factor models. Table (ref) also reports the percentage (Arb. ($\%$)) of estimated models with time-varying loadings and presenting arbitrage, namely selecting covariates violating the no-arbitrage restrictions. For that computation, both the time-invariant estimates and aOGL estimates avoid ex-ante arbitrage by construction. In line with our Monte Carlo results, the aLASSO procedure ends up with all the time-varying models estimated with arbitrage for both factor specifications. We conclude that the aOGL estimation achieves parsimony while avoiding arbitrage in time-varying factor models.

In the three first lines of Tables (ref) and (ref), we investigate the type of stock excess returns that exhibit time-variation in the their factor loadings. For both factor specifications, the longer the sample size, the more “action” is needed for the dynamics of the factor loadings. Indeed, the aOGL method selects a time-invariant model for more than 50% of the stocks excess returns exhibiting historical data smaller than 10 years, and this for both factor specifications. On the contrary, 80% of the models with the longest sample size ($\geq$ 50 years) need time-variation in their factor loadings.

In the next lines of Tables (ref) and (ref), we further look at the selected variables among the 6 in $Z_{t-1}$ and the 13 in $Z_{i,t-1}$. Across all sample sizes $T_i$, the percentages of selected variables in $Z_{t-1}$ are much higher than the percentages of selected variables in $Z_{i,t-1}$. For smaller time spans, some characteristics are never selected. Therefore, it seems that the common instruments $Z_t$ are key drivers of the time variation of the factor loadings. It is particularly true for the range $\geq$ 50y. It shows the need of including common instruments that pick up the influence of the business cycles on the factor loading dynamics in larger time spans. While there are no common instruments that are never selected, and this across all time spans, the proportions far below 100% demonstrate the need to select instruments in a data-driven way.

In Tables (ref) and (ref), we report the percentage that the 6 variables in $Z_{t-1}$ (scaled factors) and the 13 variables in $Z_{i,t-1}$ are selected through the aOGL method for both factor specifications. The percentages of selected common instruments are similar for the factors $f_{m}$, $f_{hml}$, and $f_{smb}$ shared between the two models. With the Carhart four-factor specification, we need more variables in $Z_{i,t-1}$ to describe the dynamics of the factor loadings in comparison with the Fama-French five-factor specification. In line with chaieb2020factors, the characteristics are not necessarily paired more often with their corresponding factors. The size characteristic LME is more often paired with the market factor $f_{m}$ than with the size factor $f_{smb}$. On the contrary, the momentum characteristics $r_{12,7}$ and $r_{2,1}$ are often associated with its corresponding factor $f_{mom}$. Finally, Tables (ref) and (ref) show that the conditional expectations of the factors $f_{rmw}$ and $f_{cma}$ in the Fama-French five-factor specification need less covariates than for the other factors. The variables ntis and def_spread are selected for all factors.

Let us now investigate in-sample predictability performance. As, in chaieb2020factors, we decompose the conditional expected return of asset $i$ for month $t$ for both time-varying factor specifications, as:

For such time-varying specifications, the contribution of the pricing errors $a_{i,t} - b_{i,t}^{\top}\nu_t$ is often small, revealing that the no-arbitrage restrictions are met for a vast majority of dates. When they are not, chaieb2020factors show that incorporating pricing errors, instead of only relying on $b_{i,t}^{\top} \lambda_t$ in ((ref)), helps to predict future equity excess returns. Similarly, for the time-invariant models, we decompose the unconditional expected return as:

For such time-invariant specifications, the contribution of the pricing errors $a_{i} - b_{i}^{\top}\nu$ is often large. We also consider the case of constant $\nu$ and time-varying risk premia $\lambda_t$ ($\lambda_t \& \nu$), for which we decompose the conditional expected return as

In such a hybrid model avramov2004stock, the time-variation in $E[R_{i,t}\vert {\cal F}_{t-1}]$ only comes from the time-variation in $E[f_{t}\vert {\cal F}_{t-1}]$ since $\nu$ is constant because of the no-arbitrage restrictions with constant $b_i$ and $a_i$.

To compare the prediction performance of the four estimation approaches, we compute the RMSPE of an equally-weighted portfolio for the Carhart four-factor model and Fama-French five-factor model. Equal weighting corresponds to cross-sectional averaging. chaieb2020factors also uses this weighting scheme. For that portfolio, we compute the PE by comparing the prediction made at time $t$ by each model (((ref)) and ((ref))) to the forward 12-months realized excess returns, namely the average of the realized excess returns over the next 12 months. Table (ref) reports the RMSPE, as well as the Av$(\vert \text{PE}\vert)$ and Std$(\vert\text{PE}\vert)$ for the Carhart four-factor model and Fama-French five-factor model specifications. The aOGL method performs better than its natural competitor, the aLASSO, even for that very diversified stable portfolio, where we expect differences in prediction performance to be attenuated. It is comparable in terms of the RMSPE to the $\lambda_t \& \nu$ method, with a lower Std$(\vert\text{PE}\vert)$. Figure (ref) displays the corresponding box-plots of the PE computed at each month for each method. The box-plots for the aOGL method in Figure (ref) are narrower than for the aLASSO method, and comparable for the two other methods. Those predictability improvements against the aLASSO approach provide further evidence in support for the aOGL approach advocated for the first-pass regression, so that we can incorporate model parameter restrictions to get models compatible ex-ante with the no-arbitrage restrictions. To further investigate time-varying predictability, Figures (ref) to (ref) show the forward 12-months realized excess returns for the equally-weighted portfolio and compare them with the predicted excess returns computed from ((ref)) and ((ref)) for the two methods with penalisation, respectively for the Carhart four-factor and Fama-French five-factor specifications. In both Figures (ref) and (ref), the aOGL predicted excess return paths (red plain line) overall reconcile well with the realized excess returns (black dashed line). On the contrary, the aLASSO method in Figures (ref) and (ref) does not reconcile well the predicted excess returns with the realized excess returns and sometimes predicts large negative excess returns, which is at odd with a positive reward expected from taking risks. The observed differences in the decomposition between estimates of $a_{i,t}$ (orange shaded area) and of $b_{i,t}^{\top}\mathbb{E}[f_{t}\vert\mathcal{F}_{t-1}]$ (blue shaded area) come from the selected regressors in the first pass. Since the aLASSO penalization ends up with time-varying models presenting arbitrage, we observe larger values for estimated $\hat a_{i,t}$, especially during the recession periods (gray areas) determined by the National Bureau of Economic Research (NBER). The aOGL method avoids putting covariates in estimated $\hat a_{i,t}$ that should not be there because of the no-arbitrage restrictions. Besides, the estimated path for $a_{i,t}$ is close to zero with the aOGL method as it should be if we believe that the factors are most of the time fully tradable.

Out-of-sample prediction performance

In this section, we compare the out-of-sample prediction performance for the same methods used in the previous section. Here, we compute PE but for data that never enter into model estimation. We follow a similar approach to gu2020empirical. We split the sample into two subsamples, one for training and one for testing. We estimate the models from July 1963 to December 2009 and compute PE from January 2010 to December 2019 (recent period). We repeat the same analysis for a training period from July 1963 to December 1999 and a testing period from January 2000 to December 2009 (older period). We closely follow the same setting as in the previous section, the only difference being that we separate the subsample used for estimation from the one used for prediction performance assessment.

We see that the aOGL method performs better than the aLASSO method in all cases as shown in Tables (ref) to (ref). Furthermore, the aOGL method often performs better than a time-invariant method as exhibited by the RMSPE and the lower Std$(\vert\text{PE}\vert)$. Such an advantage over time-invariant alternatives is less clear for the out-of-sample $R^2$ computed each year on the whole testing periods in Tables (ref) and (ref). On the contrary, the aOGL method keeps a strong advantage over the aLASSO method, especially for the years closer to the training periods. For both testing periods, the box-plots in Figure (ref) show that out-of-sample PE related to the portfolio excess returns for the aOGL method are located closer to zero, more symmetrically distributed, and narrower. As observed in the in-sample analysis, the aOGL method seems to perform better in terms of out-of-sample predictability as shown by the distributional behavior of the PE. We believe that the good out-of-sample performance for the portfolio comes from the diversification of the prediction errors among the single assets. We observe a similar phenomenon in forecast combinations Timmermann_2006.

Conclusions

Our empirical results show that taking explicitly into account the no-arbitrage restrictions coming from the Arbitrage Pricing Theory do help in predictive modeling of large cross-sectional equity data sets with penalisation methods. We view this approach as an example of a structural approach to big data where incorporating finance theory improves on the prediction performance of the estimated quantities. It resonates with structural approaches in panel econometrics guided by economic theory bonhomme2017keeping. In asset management and risk management, a better predictive performance of excess returns should help to better gauge time-variation in the risk-reward trade-off. In asset selection, it should help to improve performance of time-varying portfolio allocation when we use predicted excess returns as inputs. From our simulation and empirical results, we expect our procedure to perform well in out-of-sample prediction for portfolio building.

Acknowledgements

We would like to thank the Co-Editor, Associate Editor, and two referees for constructive criticism and numerous suggestions which have led to substantial improvements over the previous version. We are grateful to I.\ Chaieb, T.\ Chordia, A.-P.\ Fortin, P.\ Gagliardini, R.\ Garcia, M.\ Karemera, H.\ Langlois, A.\ Patton, S.\ Pruitt, J.\ Thimme, S.\ van den Hoff, D.\ Xiu, P.\ Zaffaroni for their helpful comments, as well as seminar participants at Queen Mary, University of Nottingham, Louvain Finance Seminar, IAE Lille, Universit\'e d'Orl\'eans, LUISS, Auburn University Statistics and Data Science Seminar, CREST-ENSAE Financial Econometrics Seminar, and participants at the $8^{th}$ days of Econometrics for Finance, 2021 SoFiE Machine Learning conference, Vienna Workshop "Econometrics of Option Markets", 2021 SFI Research days, 2021 North America Summer Meeting of the Econometric Society, $13^{th}$ Annual SoFiE Conference, 2021 IAAE Conference, 2021 EcoSta, $7^{th}$ IYFS Conference, 2021 China International Conference in Finance (CICF), 2021 EEA-ESEM, $20^\text{th}$ Workshop in Econometrics for Finance, 2021 CFE-CM Statistics, $19^\text{th}$ Paris December Finance Meeting, ML approaches Finance and Management workshop, 38$^{th}$ AFFI conference, and QFFE conference 2022. G.\ Bakalli and O.\ Scaillet were supported by the SNSF Grant $\#100018-182582$. S.\ Guerrier was supported by the SNSF Professorships Grant \#176843 and by the Innosuisse-Boomerang Grant \#37308.1 IP-ENG.