# Deep Partial Least Squares for Empirical Asset Pricing111This paper was
presented as a keynote talk at the R in Finance Conference, Chicago, 2022. The
authors would like to thank Peter Carl (Balyasny Asset Management), Francois
Cocquemas (Florida State University), and Kris Boudt (KU Leuven), for useful
discussions and input on this paper. 222The source code and data for this
paper is available at https://github.com/Kem975/Deep_Partial_Least_Squares.

Matthew Dixon333Email<EMAIL_ADDRESS>
Department of Applied Mathematics
Illinois Institute of Technology
Nicholas G. Polson
Booth School of Business
University of Chicago
Kemen Goicoechea
Department of Applied Mathematics
Illinois Institute of Technology

###### Abstract

We use deep partial least squares (DPLS) to estimate an asset pricing model
for individual stock returns that exploits conditioning information in a
flexible and dynamic way while attributing excess returns to a small set of
statistical risk factors. The novel contribution is to resolve the non-linear
factor structure, thus advancing the current paradigm of deep learning in
empirical asset pricing which uses linear stochastic discount factors under an
assumption of Gaussian asset returns and factors. This non-linear factor
structure is extracted by using projected least squares to jointly project
firm characteristics and asset returns on to a subspace of latent factors and
using deep learning to learn the non-linear map from the factor loadings to
the asset returns. The result of capturing this non-linear risk factor
structure is to characterize anomalies in asset returns by both linear risk
factor exposure and interaction effects. Thus the well known ability of deep
learning to capture outliers, shed lights on the role of convexity and higher
order terms in the latent factor structure on the factor risk premia. On the
empirical side, we implement our DPLS factor models and exhibit superior
performance to LASSO and plain vanilla deep learning models. Furthermore, our
network training times are significantly reduced due to the more parsimonious
architecture of DPLS. Specifically, using 3290 assets in the Russell 1000
index over a period of December 1989 to January 2018, we assess our DPLS
factor model and generate information ratios that are approximately 1.2x
greater than deep learning. DPLS explains variation and pricing errors and
identifies the most prominent latent factors and firm characteristics.

Key Words: Deep Learning, Partial Least Squares, Shrinkage, Empirical Asset
Pricing, Conditional Latent Factor Models, Non-Linear Risk.

##  1 Introduction

In this paper, we develop a deep factor methodology known as deep partial
least squares (DPLS) as a projection based dimensionality reduction technique
for empirical asset pricing regressions which, combined with deep learning,
provides improved predictive performance on financial panel data while also
providing information on the risk factor premia. Our DPLS approach is shown to
be a natural, composable, and powerful approach to projecting both the
predictors and the responses onto a relatively small set of orthonormal latent
variables so as to maximize the covariance between the projected predictors
and the projected responses. In contrast to principal component dimension
reduction, PLS & DPLS projections are dependent on the response and thus
leverage a wider dataset of firm charactersistcs for predicting returns.
Polson et al. (2021) find evidence that DPLS shows superior predictive
performance on non-financial cross-sectional data. Dixon and Polson (2020)
show that deep learning identifies the most important interaction effects and
factors. DPLS builds on this by boosting the covariance between the response
and predictor, while simultaneously using the same statistical loadings
generated by PLS. Furthermore, because DPLS is based on deep learning, it is
able to automatically capture interaction effects.

Partial Least Squares (PLS) is a wide class of methods for modeling relations
between sets of observed variables by means of latent variables. It is
comprised of regression and classification tasks as well as dimension
reduction techniques and modeling tools. The underlying assumption of all PLS
methods is that the observed data is generated by a system or process which is
driven by a smaller number of statistical latent variables.

Across hedge funds, asset management, and proprietary trading firms, it is
commonplace to use supervised learning for macroeconomic and financial
forecasting by performing ”feature engineering” — expanding the original set
of covariates to include terms such as interactions, dummy variables or
nonlinear functionals of the original predictors (Hoadley, 2001; Dobrev and
Schaumburg, 2013; Dixon et al., 2017). It’s also increasingly common to
aggregate multiple data sources for financial forecasting, combining
traditional financial data such as stock prices with novel signals such as
news sentiment (Dierckx et al., 2020; Algaba et al., 2021). Our goal is to
employ supervised dimension reduction techniques which, when combined with
deep learning, not only leads to a more parsimonious architecture with
improved out-of-sample performance, but enable factor models to be properly
seated in an equilibrium asset pricing framework in which the statistical
factors serve as latent risk factors and a non-linear risk factor structure is
represented by the network architecture. Our approach corresponds to the use
of a non-linear stochastic discount factor (see Cochrane (2011)).

The rest of the paper is outlined as follows. Section 1.1 outlines the general
connections with previous work and highlights the contribution of this paper.
Section 2 surveys the related works on factor modeling works in more detail
culminating in the introduction of our non-linear factor model, best viewed as
a non-linear generalization of the BARRA factor model. Section 3 describes our
general modeling methodology using deep learning within partial least squares
for multi-index regression. Following Polson et al. (2021), we provide a
nonlinear extension of partial least squares by using deep layers to model the
relationship between the scores generated in PLS. Then in Section 3.4, we
present the DPLS factor model and some of its variants. Section 4 demonstrates
the application of our DPLS factor model to equity data consisting of monthly
returns and firm characteristics. Finally, Section 5 concludes with directions
for future research.

###  1.1 General Connection to Previous Work

Within the investment management industry, asset managers seek novel
predictive firm characteristics to explain anomalies which are not captured by
classical capital asset pricing and factor models. Recently a number of
independent empirical studies, have shown the importance of using a higher
number of, often highly correlated, predictors related to firm characteristics
and other common factors (Moritz and Zimmermann, 2016; Harvey et al., 2015; Gu
et al., 2018; Feng et al., 2018). In particular, Gu et al. (2018) analyze a
dataset of more than 30,000 individual stocks over a 60 year period from 1957
to 2016, and determine over 900 correlated baseline signals.

Despite being widely used in the asset management industry, it is becoming
increasingly well documented that ordinary linear regression is inadequate for
high dimensional factor models. The frequently encountered problem of nearly
collinear regressors can be addressed using regularization to provide
shrinkage estimation (Hahn et al., 2013; Dobrev and Schaumburg, 2013; Cheng et
al., 2016). Dobrev and Schaumburg (2013) develops a Regularized Reduced Rank
Regression model, a special case of which generalizes principal component
regression by applying reduced rank rather than linear regression to the
principal components of the regressors, thereby disentangling the forecasting
factors driving the outcomes from the factor structure in the predictors.

Moritz and Zimmermann (2016); Gu et al. (2018); Chen et al. (2019) highlight
the inadequacies of OLS regression in variable selection over high dimensional
datasets — in particular the inability to capture outliers. In contrast, deep
neural network can explain more structure in stock returns because of their
ability to fit flexible functional forms with many covariates, without apriori
knowledge of the interaction effects. Feng et al. (2018) demonstrates the
ability of a three-layer deep neural network to effectively predict asset
returns from fundamental factors. Gu et al. (2021) apply auto-encoders to fit
a latent factor conditional asset pricing model and demonstrate that the
approach not only out-performs other factor models but is suitable for no-
arbitrage pricing. However, techniques which project on statistical factors
are difficult to economically interpret and have limited hedging utility.

Dixon and Polson (2020) demonstrate a three-layer deep neural network approach
for fundamental factor modeling which not only outperforms OLS regression and
LASSO but also provides interpretability –specifically, the importance of the
predictors and interactions can be ranked and the uncertainty in predictor
sensitivities can be estimated with bootstrap sampling techniques. Dixon and
Polson (2022) extend their methodology to a Bayesian estimation framework to
provide uncertainty quantification and a predictive distribution of excess
returns.

There are a number of advantages of DPLS methodology over traditional deep
learning. First, PLS uses linear methods to extract predictors and weight
matrices. Due to a fundamental estimation result of Brillinger (2012) we show
that these estimates can provide a stacking block for our DPLS model.
Moreover, traditional diagnostics (bi-plots, scree-plots) to diagnose the
depth and nature of the activation functions allow us to find a parsimonious
model architecture. Frank and Friedman (1993) provides a Bayesian shrinkage
interpretation of PLS. This permits the discrimination between unsupervised
and supervised learning method. In other words, based on the value of
shrinkage factors, one can deduce whether a supervised learning method such as
DPLS is necessary or whether an unsupervised dimension reduction method such
as PCA is adequate.

One of the remarkable discoveries by Polson et al. (2021) in fitting DPLS
regression models is that the PLS projection can be entirely separated from
the training of the deep network, thus enabling standard deep neural networks
algorithms, such as stochastic gradient descent, to be used without
modification. Drawing on a seminal result by Brillinger (1977) which paved the
way for the landmark paper by Brillinger (2012) and later Naik and Tsai
(2000), we formalize the critical observation that under some limiting
assumptions, the projection directions in PLS are invariant to the functional
form of the regression between score matrices, up to a constant of
proportionality.

Our theoretical development parallels Brillinger (2012), namely that PLS
regression is a strongly consistent estimator to the regression coefficient
even when there exists a unknown non-linear relationship between the
predictors and the responses. We invert this line of reasoning and instead
show that the linear PLS estimate is strongly consistent with the score
regression coefficients when the regression on scores is a (known) non-linear
function. This unlocks PLS regression as a composable method – non-linear
functional approximations can be composed with the PLS projection matrices
without violating the consistency result. More precisely, one can use still
use the same projections and x-loadings estimated by PLS, under a linear
relationship between the input and output scores, while subsequently using
deep learning to improve the predictive performance of the y-scores.

In the context of empirical asset pricing, the DPLS method is shown to deliver
a key point of departure from Gu et al. (2018); Chen et al. (2019); Kelly et
al. (2019) in enabling the modeling of a _non-linear factor structure_ which
corresponds to an assumed _non-linear stochastic discount factor_.
Furthermore, we develop interpretative techniques for DPLS which show both
firm characteristic and latent factor sensitivities and interaction effects
are provided by the DPLS factor model.

As a further contribution, we show how our DPLS factor models not only provide
a non-linear generalization of BARRA fundamental factor models but also
briefly discuss how they provide a non-linear, non-stationary extension of
dynamic factor models (Stock and Watson, 2010).

On the empirical side, we develop a no-arbitrage DPLS factor model with 50
factors for 3290 Russell 1000 indexed stocks over an approximate 30 year
period and compare performance with LASSO and pure deep learning based
fundamental factor models. Furthermore, we assess the risk-reward trade-off of
our factor model and identify the implications of including a non-linear
factor structure.

##  2 Factor Modeling

Our goal is to explain the differences in the cross-section of returns for
individual stocks using a latent factor no-arbitrage asset pricing model,
whereby any expected gain in asset returns over the risk-free rate is solely
attributed to factor risk exposure. But rather than this returns model being
linear in factors, or equivalently, assuming linear dependence of a stochastic
discount factor on the factors, we seek to introduce non-linearity. The
rationale for non-linearity in factor modeling is well motivated in the
literature and can be reasoned from Stein’s lemma, that the non-linearity in
the dependency of the returns on the factors arises precisely when the asset
returns and factors are non-Gaussian.

While the practice of applying non-linear parametric functions, such as deep
learning, to factor models is conceptually straightforward, the challenge
becomes in how to go beyond a mere mechanical exercise in learning a non-
linear function of predictor variables that maximizes the out- of-sample
explanatory power for realized excess returns of stocks, to characterizing the
risk-return tradeoff of assets built from our factor model, under the
equilibrium asset pricing principle that characteristic-based predictability
are determined through risk exposures to factors alone (Kelly et al., 2019).
Thus, in informal terms, in order for factor models to be correctly seated in
the framework of equilibrium asset pricing, it must resolve the predictability
attributed to the factor exposures and not conflate this with mis-pricing.
Furthermore, they must be no-arbitrage factor models and hence provide no
intercept when regressed against excess returns.

The need for both explanatory power from a set of predictor variables, taken
to be lagged firm characteristics, and risk-return properties inevitably leads
to model design trade-offs. For example some machine learning approaches to
factor modeling yield strong predictive power but yield no insight into risk
factor attribution and therefore have questionable economic interpretability
and practical utility. On the other hand, other approaches deploy machine
learning in such a way as to guarantee that the characteristic-based
predictability is attained solely through factor risk exposures while
accounting for non-linearity between the factor loadings and the firm
characteristics. However, those approaches that do resolve the risk factors,
still rely on a linear relationship between the returns and factor loadings,
or equivalently the stochastic discount factor remains linear in the factors
(Chen et al., 2019).

Then there are questions of how best to reduce the set of predictor variables
to a smaller number of latent factors to manage the risk-reward tradeoff. Such
factors should ideally carry economic interpretability and be representable by
factor mimicking portfolios. We can broadly separate the main approaches to
non-linear factor models into the following three categories:

  1. 1.

Forecasting models which predict the excess returns directly as a non-linear
function of the lagged firm characteristics. This essentially amounts to a
pure exercise in prediction of excess returns given the lagged firm
characteristics, without assessing the risk/return tradeoff in terms of the
latent factors (see Gu et al. (2018); Dixon and Polson (2020)).

  2. 2.

Linear Conditional Latent Factors preserve a linear relationship between the
excess returns and the factors – amounting to the use of a linear stochastic
discount factor in the factors – but introducing the non-linearity between the
conditional factor loadings and the lagged firm characteristics (see, for
example, Kelly et al. (2019); Chen et al. (2019); Gu et al. (2021)). These
approaches allow for latent factors and time-varying loadings which “ _can
explain individual stock returns and estimate the best performing factors
without taking an a prior stand on what the factors are_ ” Kelly et al.
(2019). This category of approaches is thus suitable for asset pricing as it
provides insight into the factor risk premia by conditioning on firm
characteristics.

  3. 3.

Non-linear Conditional Latent Factors model the excess returns directly as a
non-linear function of the lagged firm characteristics, but allow for the
estimation of the conditional factor loadings and the latent factors under a
non-linear stochastic discount factor in the factors. Hence, the approach is
suitable for asset pricing and it operates under the less restrictive
assumption that the stochastic discount factor is linear in the factors.
Equivalently, as per Stein’s lemma, such an approach does not assume that the
asset returns and factors are Gaussian.

Our DPLS factor model bridges these first two categories to arrive at a model
in the final category in its most general form. We stop short of claiming that
the modeling approach recovers the second category of approaches, as there are
nuanced issues in how the non-linearity is represented, the estimation of the
latent factors, and the handling of the intercepts. Specifically, our approach
is to use deep projected least squares regression to identify _non-linear
factor structure_ , that is the asset returns are attributed to non-linear
risk factors, with non-linear conditional factor loadings in the lagged firm
characteristics, rather than linear risk factors under non-linear conditional
factor loadings. Moreover, the approach is suitable for no-arbitrage pricing
as the intercept captures the portion of the model returns not attributed to
the risk factors. We thus build on several ideas which culminate from the
extensive factor model literature, some of the most closely related works we
now detail here.

Gu et al. (2018) present a static predictive model for the excess returns at
time $t$, $\mathbf{r}^{e}_{t}:=[\mathbf{R}^{e}]_{t}$, which learns a flexible
parameterized function, $h_{\theta}(\cdot)$ of a set of predictors – lagged
firm characteristics, $\mathbf{z}_{t}$, observed at time $t$:

$r^{e}_{i,t+1}=h_{\theta}(z_{i,t})+\epsilon_{i,t+1},\leavevmode\nobreak\
i=1,\dots,N_{t},\leavevmode\nobreak\ t=1,\dots,T,$

where $\bm{\epsilon}_{t}:=[E]_{t}$ are the idiosyncratic zero-mean Gaussian
errors. Dixon and Polson (2020) present a time dependent parameterization of
this model which is trained period-by-period. The advantage of the time
dependent model is that it needn’t assume stationarity in the returns and firm
characteristics. Rather, at each period, a new model is fitted, which only
relies on data from the current period and not over a questionably relevant
history. It does however result in fewer training observations to fit the
model and hence shrinkage, or other forms of model reduction, are needed to
prevent over-fitting of highly parameterized models such as deep neural
networks. Equivalently, Bayesian inference under a choice of prior on the
parameters can reduce the number of training samples needed (see Dixon and
Polson (2022) for further details of a Bayesian deep fundamental factor
model).

These predictive models, however, offer no insight into the latent factor
structure and therefore offer no insight into whether the
characteristics/expected return relationship is driven by compensation for
exposure to latent risk factors (Kelly et al., 2019). Specifically, there
appears to be mechanism for identifying how the characteristics proxy for
loadings on common risk factors.

Gu et al. (2021) present a conditional latent factor model for excess asset
returns, in which the conditional and dynamic $K-$factor loadings
$\bm{\beta}(z_{i,t})$ are a non-linear parametric function of $p$ lagged firm
characteristics, $z_{i,t}$. Additionally the excess asset returns are linear
in the $K$-latent factors, $\mathbf{f}_{t}$, which are in turn a linear
combination of all excess returns in the universe. The model is fitted period-
by-period, i.e. $t=1,\dots,T$ as follows:

$\displaystyle r^{e}_{i,t+1}$ $\displaystyle=$
$\displaystyle\bm{\beta}(z_{i,t})^{\top}\mathbf{f}_{t+1}+\epsilon_{i,t+1},\leavevmode\nobreak\
i=1,\dots,N_{t},$ (1) $\displaystyle\bm{\beta}(z_{i,t})$ $\displaystyle:=$
$\displaystyle g_{\theta_{t}}(z_{i,t}),$ (2) $\displaystyle f_{k,t}$
$\displaystyle:=$
$\displaystyle\mathbf{w}_{k,t}^{\top}\mathbf{r}^{e}_{t},\leavevmode\nobreak\
k=1\dots K.$ (3)

The advantage of their approach is three-fold: (i) asset return predictability
derives entirely from the lagged firm characteristics via the exposures to the
risk factors; (ii) the factors aren’t pre-determined heuristically by a two-
step procedure, as in the Fama-French model, but rather the factor structure
is determined statistically in a single-step, taking a data driven approach to
factor structure derivation and not relying on heuristics for characterizing
the factor-structure; and (iii) the factors can be interpreted from how the
dynamical loadings map to firm characteristics. For example, in any given
period, their model can characterise which firm characteristics are the most
prominent across all assets in any factor loading.

A technical development in Gu et al. (2021) is the presentation of an auto-
encoder like neural network architecture which is customized to represent the
structure of the model. Specifically a feed-forward architecture for the firm
characteristics yields a hidden layer output representing the factor loadings
and this is subsequently multiplied by latent factors which are given by a
linear feedforward layer, with $N_{t}$ excess asset returns as inputs. While
this network exhibits the signature bottleneck of an autoencoder, it is
technically not an auto-encoder. This is because, like PCA, autoencoders are
unsupervised learning methods which rely solely on returns to identify the
latent factor structure. They do not leverage conditioning variables to
identify this factor structure.

However, in seeking to modify this architecture to include firm characteristic
data, they have essentially arrived at what partial least squares regression
was exactly designed for, namely supervised dimension reduction by accounting
for the returns and leveraging the lagged firm characteristics to identify the
latent factor structure. Hence, although Gu et al. (2021) serve as the
inspiration for our work, we argue that the deep partial least squares
regression architecture is the natural choice for non-linear conditional
latent factor models.

Another important distinguishing aspect of factor models is their approach to
variable selection and regularization. The literature seems to be divided in
statistical variable selection and regularization approaches versus the use of
financial theory or “model based” variable selection and regularization. A
prime example of statistical approaches would be LASSO penalized linear models
which use shrinkage and variable selection to manage high dimensionality by
forcing the coefficients on most regressors near or exactly to zero. This can,
however, produce sub-optimal forecasts when predictors are highly correlated.
A simple example of this problem is a case in which all of the predictors are
equal to the forecast target plus an iid noise term. In this situation,
choosing a subset of predictors via LASSO penalty is inferior to taking a
simple average of the predictors and using this as the sole predictor in a
univariate regression. However, other approaches to variable selection are
commonly used in factor modeling. PCA, PCR, and PLS are all examples of
projection methods which rely on a shrinkage mechanism which shrinks away from
the origin for certain eigen-directions. The resulting shrinkage factors along
each of the eigen-directions, referred to as principle components in PCA and
PCR, or scale factors in PLS, are chosen to maximize the explained covariance
either in the predictors or between the predictors and the responses
respectively. These static approaches appear widely in the literature, most
notably the asymptotic PCA (APCA) method (Connor and Korajczyk, 1988) which is
suitable when $N>T$, and the use of PCA (Swanson, 2014) when $N<T$. As a panel
regression, $\mathbf{R}^{e}=BF+E$, Stock and Watson (2010) and Bai and Ng
(2002) apply PCA to the covariance of asset returns:

$\bar{\mathbf{R}}^{e}=\hat{P}\Lambda\hat{Q}^{\top}+\hat{E},$

where $\hat{P}$ and $\hat{Q}$ are the left and right singular vectors of the
covariance matrix, representing the factor loadings, and $\Lambda$ is a
diagonal matrix with eigenvalues, sorted in descending order, representing the
factors (i.e. principle components). The first latent factor corresponds to
the direction that captures the maximum variance in the asset returns, the
second factor corresponds to the direction of maximum variance under a
complementary projection in the remaining orthogonal subspace, and so forth.
This parsimony of factors is appealing as it can explain asset return
covariance by a small number of systematic factors across the market.

These approaches suffer from many shortcomings – aside from violating the no-
arbitrage principle unless the PCA methodology is modified to use uncentered
second moments of excess returns– the main shortfalls being that the latent
factors poorly explain the outliers in asset returns. As put succinctly by
Kozak et al. (2018), “ _It is perfectly possible for a [principle component]
factor to be important in explaining return variance but to play no role in
pricing_ ”. Another point of distinction is that the factors are purely
statistical and not macro-financial, hence such risks are difficult to
interpret in a macro-financial model.

A key development relevant to this paper is Kelly and Pruitt (2015) , who
apply PLS as a three-pass regression filter to successfully filter out the
impact of irrelevant risk factors, a problem that the method of principal
components alone cannot be guaranteed to do. However, unlike their three-step
regression, our approach is a period-by-period single step model and, like Gu
et al. (2021), our betas are dynamic. So while our latent risk factors may
vary period-to-period, we avoid data stationarity assumptions which are needed
by Kelly and Pruitt (2015).

Another line of reasoning is pursued by Lettau and Pelger (2020), who propose
a new method for estimating latent asset pricing factors that fit the time-
series and cross-section of expected returns. Their estimator generalizes PCA
by including a penalty on the pricing error in expected return leading to
statistical risk factors which depend only on small set of firm
characteristics. Their work sheds light on the “factor zoo” problem, namely
which risk factors are really important and which factors are subsumed by
others (Cochrane, 2011). While our model is period-by-period cross-sectional
and does not mirror the experimental design of Lettau and Pelger (2020), we do
seek to identify how, in any period, out-of-sample predictions of asset
returns are attributed to the statistical risk factors and which firm
characteristics, combined with their interaction effects, are most prominent.

Another modernized approach, which generalizes PCA to allow both dynamic
loading and factor structure derived from firm characteristics, is the
instrumented PCA (IPCA) factor model Kelly et al. (2019), which for ease of
exposition, we show in its restricted form:

$\displaystyle r^{e}_{i,t+1}$ $\displaystyle=$ $\displaystyle
z_{i,t}^{\prime}\Gamma f_{t+1}+\epsilon_{i,t+1},\leavevmode\nobreak\
i=1,\dots,N_{t},\leavevmode\nobreak\ t=1,\dots,T,$ (4)

where $\Gamma\in\mathbb{R}^{p\times K}$ projects the lagged firm
characteristics on to the $K<p$ latent factors and is static. Note that the
loadings are now dynamic due to the reliance on the lagged firm
characteristics. This model allows for a potentially large number of
characteristics to be mapped on to a smaller set of common latent risk factors
that explain anomalies in asset returns. The projected characteristics, or
more precisely, a linear combination of characteristics, thus proxy for
loadings on each of the common risk factors and admit the intuitive
interpretation that higher loadings result in higher expected returns to
compensate for higher exposures, i.e. increased factor risk premia. More
importantly, the most important characteristics associated with each risk
factor can be estimated from observed excess returns and firm characteristics
without, apriori, the use of heuristic mappings between them.

Our experimental design broadly follows the cross-sectional modeling approach
used by (Fama and MacBeth, 1973) and BARRA factor models (see Rosenberg and
Marathe (1976); Carvalho et al. (2012)) are appealing because of their
simplicity and their economic interpretability, generating tradable
portfolios. The BARRA model is able to capture linear effects of macroeconomic
events on individual securities. The choice of factors are microeconomic
characteristics - essentially common factors, such as industry membership,
financial structure, or growth orientation (Nielsen and Bender, 2010).

Under the assumption of homoscedasticity, factor realizations can be estimated
in the BARRA model by ordinary least squares regression444The BARRA factor
model is often presented in the more general form with heteroscedastic error
but we do not consider the non-linear extension here.. OLS regression exhibits
poor expressability and relies on the Gaussian errors being independent of the
regressors. Generalizing to non-linearities and incorporating interaction
effects is a harder task.

The BARRA fundamental factor model555Specifically, the predictive form of the
BARRA model is conventionally used for risk modeling, however, we shall use
adopt a non-linear version of this model for a stock selection strategy based
on predicted returns. is a single-index cross-sectional model which assumes
that the factor loadings are proxied by the firm characteristics. The model
expresses the linear relationship between $K$ fundamental factors and the
asset’s returns:

$\mathbf{r}_{t+1}=B_{t+1}\mathbf{f}_{t+1}+\mathbf{\epsilon}_{t+1},\leavevmode\nobreak\
t=1,\dots,T,$ (5)

where $B_{t}=[\mathbf{1}\leavevmode\nobreak\ |\leavevmode\nobreak\
{\bm{\beta}}_{1}(t)\leavevmode\nobreak\ |\cdots|\leavevmode\nobreak\
{\bm{\beta}}_{K}(t)]=[\mathbf{1},Z_{t-1}]$ is the $N\times K+1$ matrix of
known factor loadings (betas):
$\beta_{i,k}(t):=\left(\bm{\beta}_{k}\right)_{i}(t)$ is the exposure of asset
$i$ to factor $k$ at time $t$ and assumed to proxy the $k^{th}$ firm
characteristic. The loadings therefore correspond to asset specific attributes
such as market capitalization, industry classification, style classification.
$\mathbf{f}_{t}=[\alpha_{t},f_{1,t},\dots,f_{K,t}]$ is the $K+1$ vector of
unobserved factor realizations at time $t$, including $\alpha_{t}$.
$\mathbf{r}_{t}=\frac{P_{t}}{P_{t-1}}-1$ is the $N-$vector of single-period
asset returns observed at time $t$. The errors are assumed independent of the
factor realizations $\rho(f_{i,t},\epsilon_{j,t})=0,\forall i,j,t$ with
heteroschedastic Gaussian error,
$\mathbb{E}[\epsilon^{2}_{j,t}]=\sigma^{2}_{j}$. Then, in principle, factor
realizations can be estimated with OLS regression. Note that there are more
general estimation frameworks for estimating factors under heteroschedastic
error (see e.g. Korteweg and Polson (2009); Tsay (2010)). The BARRA model can
be viewed as a special case of IPCA when the projection is the identity matrix
(Kelly et al., 2019).

Following Dixon and Polson (2020), we can generalize the BARRA model by
considering a period-by-period non-linear cross-sectional fundamental factor
model of the form

$\mathbf{r}_{t+1}=F_{t+1}(Z_{t})+\mathbf{\epsilon}_{t+1},\leavevmode\nobreak\
t=1,\dots,T,$ (6)

where $\mathbf{R}_{t+1}$ are single-period asset returns,
$F_{t+1}:\mathbf{R}^{K}\rightarrow\mathbf{R}$ is a differentiable non-linear
function that maps the $i^{th}$ row of $Z_{t}$ to the $i^{th}$ asset return at
time $t+1$. The map is assumed to incorporate a bias term so that
$F_{t+1}(\mathbf{0})=\bf{\alpha}_{t+1}$. In the special case when
$F_{t+1}(Z_{t})$ is linear, the map is $F_{t+1}(Z_{t})=Z_{t}\mathbf{f}_{t+1}$.
$\mathbf{\epsilon}_{t+1}$ are zero mean i.i.d. error variates.

The forecasted returns, $\hat{r}_{t+1}$ are given by the conditional
expectation, $\hat{r}_{t+1}=\mathbb{E}[r_{t+1}\leavevmode\nobreak\
|\leavevmode\nobreak\ Z_{t}]=F_{t+1}(Z_{t})$. Hence, we seek to predict
returns $\mathbf{r}_{t+1}$ over the horizon $[t,t+1]$ from the lagged firm
charactertistcs $Z_{t}$, which are by construction measurable at time $t$.
This model can hence be mapped on to a growing body of research on returns
prediction using neural networks. One appealing aspect of this approach is
that stationarity of the factor realizations is not required since we only
predict one period ahead and, in contrast to time series models, then retrain
in each subsequent period using the latest observations only.

In Section 3.4, we show how DPLS is applied to this model to project the firm
characteristics onto a smaller set of latent risk factors with resulting non-
linear factor structure in the learned model and characterization of the risk
factor premia.

###  2.1 Stochastic discount factors

For completeness, we comment on our modeling approach from an entirely
different line of reasoning rooted in arbitrage pricing theory. Cochrane
(2011) shows how the one-factor CAPM model of Sharpe can be obtained as a
result of consumption utility maximization with the Markowitz quadratic
utility function, and under the assumption that the market is at equilibrium
and all investors are identical (and the same as the representative investor)
and hold the market portfolio. Furthermore, Cochrane (2011) also shows how
different asset pricing models can be equivalently interpreted as particular
models for the stochastic discount factor (SDF) (an ‘index of bad times’)
$m({\bf f}_{t},t)$ where ${\bf f}_{t}$ is a set of observed or inferred
dynamic variables such as market and sector returns, macroeconomic variables
etc. The SDF $m({\bf f}_{t},t)$ does _not_ depend on characteristics of
individual assets such as stocks, because with this theory the SDF is
_universal_ for all assets including stocks, bonds, derivatives etc.

|
---|---
(a) Excess returns | (b) Book-to-Price
|
(c) Total asset growth | (d) Cash flow volatility
Figure 1: Cross sectional excess returns and fundamental factors across stocks
in the Russell 1000 index in February 2013. (a) The monthly excess asset
returns (Skewness 0.96, Kurtosis 5.37); (b) Book-to-price ratios (Skewness
0.71, Kurtosis 0.58); (c) Total asset growths (Skewness -3.21, Kurtosis 40.28)
(d) Cash flow volatilities (skewness 7.76, Kurtosis 86.37).

While less visited than the above literature on estimating risk premia, there
is yet another line of reasoning that focuses on estimation of risk prices,
i.e., the extent to which the factor associated with a characteristic helps
price assets by contributing to variation (Kozak et al., 2018; Chen et al.,
2019). The former authors arrive at a characteristics-sparse stochastic
discount factor (SDF) representation that is linear in only a few such factors
and conclude that there simply isn’t enough redundancy across the factors to
justify only a handful of factors, e.g. FF3, FF4, FF5 etc. (Chen et al.,
2019)’s approach to enforcing no-arbitrage factor modeling is predicated on
the notion that the mere existence of a SDF $m({\bf f}_{t},t)$ itself is
critically dependent on the absence of arbitrage, once this assumption is
accepted, any model building amounts to designing a parameterized model
$m_{\theta}({\bf f}_{t},t)$ where $\theta$ is a set of trainable parameters,
and training the model on available data. Chen et al. (2019) introduced
Generative Adversarial Networks (GANs) to incorporate the no-arbitrage
condition as part of the neural network algorithm and estimate a linear
stochastic discount factor that explains all stock returns from the
conditional moment constraints implied by no-arbitrage. For ease of
exposition, here we drop the $t$ subscript on the SDF, returns, and factors.

It is well known, as a consequence of Stein’s lemma, that any nonlinear model
of the stochastic discount factor, $m=g(\mathbf{f})$, with a Lipschitz
continuous $f$, can be turned into a linear model
$m=a+\mathbf{b}^{\top}\mathbf{f}$ by assuming normal returns, $R$, and normal
factors, $\mathbf{f}$ (see pg. 154 of Cochrane (2011)). Such a linear
stochastic discount factor corresponds to a linear factor beta pricing model
$\mathbb{E}[R]=\alpha+\bm{\beta}^{\top}\bm{\lambda}$, where
$\bm{\lambda}=\mathbb{E}[\mathbf{f}]$.

However, as demonstrated by the cross-sections of all Russell index assets on
February 2013 in Figure 1, the assumption of normality on monthly excess
returns and factors is too strong and we therefore seek a non-linear
stochastic discount factor corresponding to a non-linear factor beta pricing
model, $h(\mathbf{\bm{\lambda}})$. We show that expressing the stochastic
discount factor $m$ in terms of a non-linear function, $g(\mathbf{f})$, in the
factors $f$ is equivalent to the non-linear factor beta pricing model.

######  Theorem 1 (Non-linear discount factors).

Given the non-linear stochastic discount factor

$m=g(\mathbf{f}),\leavevmode\nobreak\ g:\mathbb{R}^{K}\rightarrow\mathbb{R}\in
C^{\infty}(\mathbb{R}^{p}),\leavevmode\nobreak\ 1=\mathbb{E}[m]$

$\exists$ a function $h:\mathbb{R}^{K}\rightarrow\mathbb{R}\in
C^{\infty}(\mathbb{R}^{p})$ s.t. $\mathbb{E}[R]=h(\mathbf{\bm{\lambda}}).$

See Appendix C for the proof.

######  Remark 1.

Conversely, given a $h(\mathbf{\bm{\lambda}})$, $\exists$ a $g(\mathbf{f})$.

######  Remark 2.

Note for, avoidance of doubt, that we do not estimate the SDF in this paper as
in Kozak et al. (2018); Chen et al. (2019). Rather we reason, on the basis of
this theorem, that if we estimate risk premia via a non-linear function,
$h(\mathbf{\bm{\lambda}})$, then there is a corresponding non-linear SDF. The
economic implications of a corresponding utility function being non-Markowitz
quadratic are beyond the scope of this paper. However, we comment in passing
that many constrained utility functions used by practitioners to incorporate
exposure limits, maximum drawdown limits, shorting selling limits, CVaR
penalized expected returns etc, depart from the mean-variance utility function
and motivate a more flexible form of the SDF than linear in the discount
factors.

##  3 DPLS Regression

Informally, our goal is to jointly project the excess asset returns and firm
characteristics with projection matrices $Q_{t}$ and $P_{t}$ on to a pair of
components $(U_{t},V_{t})$ respectively such that:

$\mathbf{r}^{e}_{t}=U_{t}Q_{t}+\mathbf{\epsilon}_{t}^{r},\leavevmode\nobreak\
Z_{t-1}=V_{t}P_{t}+\mathbf{\epsilon}_{t}^{z},$

in a such a way that the prediction of excess asset returns

$\hat{\mathbf{r}}^{e}_{t+1}=\mathbb{E}[\mathbf{r}^{e}_{t+1}\leavevmode\nobreak\
|\leavevmode\nobreak\ Z_{t}]=\hat{U}_{t}Q_{t}$

is based on using the learned map $\hat{U}_{t}=g_{\theta}(V_{t})$ in the
projected space. The approach is thus using conditioning information from the
lagged firm characteristics and the excess asset returns to arrive at a non-
linear representation of the latent risk factors.

To formalize this method, we begin by reviewing linear PLS regression in the
more general setup before introducing the DPLS method of Polson et al. (2021).
Then in Section 3.4 we introduce a deep PLS factor model, where the factor
loadings, $V_{t}$, are interpreted as risk premia and the latent risk factors
are given by the gradient of the fitted function, $g_{\theta}(V_{t})$, in the
projected space.

In general, PLS regression uses both $X\in\mathbb{R}^{N\times p}$ and
$Y^{N\times q}$ to calculate the projection, further it simultaneously finds
projections for both input $X$ and output $Y$, making it applicable to more
general problems with multi-index output vectors as well as high dimensional
input vectors. PLS searches for a set of latent vectors (i.e. components) that
perform a simultaneous decomposition of $X$ and $Y$ with the constraint that
these components explain as much as possible of the sampled covariance between
$X$ and $Y$.

PLS hence finds projection directions which maximize the sample covariance
between $X$ and $Y$— the resulting projections $U$ and $V$ for $Y$ and $X$,
respectively are called score matrices, and the projection matrices $P$ and
$Q$ are called loadings.

Figure 2: An overview of the PLS (DPLS) method. PLS uses a linear
approximation between U and V, whereas DPLS uses a deep neural network
approximation (red line in regression plot).

The input-score $V$, has $K$ columns, one for each “feature” and $K$ is chosen
via cross-validation. The key principle is that $V$ is a good predictor of
$U$, the output-scores. These relations are illustrated in Figure 2 and
summarized by the equations below for the general multivariate PLS model:

$\displaystyle Y$ $\displaystyle=UQ+E^{y},$ $\displaystyle X$
$\displaystyle=VP+E^{x},V^{\top}V=I,$

Here $Q\in\mathbb{R}^{K\times q}$ and $P\in\mathbb{R}^{K\times p}$ are low-
rank projection (loading) matrices and $V,U\in\mathbb{R}^{N\times K}$ are
projections of $X$ and $Y$ respectively. The columns of $V$ are orthonormal
and referred to as the “latent vectors”.

DPLS uses a feedforward neural network to regress the output-score on the
input-scores:

$U=G_{W,b}(V)+E,$

where $G_{W,b}$ is a deep neural network with $L$ layers, that is, a super-
position of univariate semi-affine functions,
$\sigma^{(\ell)}_{W^{(\ell)},b^{(\ell)}}$, to give

$U:=G_{W,b}(V)=(\sigma^{(L)}_{W^{(L)},b^{(L)}}\circ\dots\circ\sigma^{(1)}_{W^{(1)},b^{(1)}})(V),$

and the unknown parameters are a set of weight matrices
$W=(W^{(1)},\dots,W^{(L)})$ and a set of bias vectors
$b=(b^{(1)},\dots,b^{(L)})$. Any weight matrix
$W^{(\ell)}\in\mathbb{R}^{m\times n}$, can be expressed as $n$ column
m-vectors
$W^{(\ell)}=[\mathbf{w}^{(\ell)}_{,1},\dots,\mathbf{w}^{(\ell)}_{,n}]$. We
denote each weight as $w^{(\ell)}_{ij}:=\left[W^{(\ell)}\right]_{ij}$.

The $\ell^{th}$ semi-affine function is itself defined as the composition of
the activation function, $\sigma^{(\ell)}(\cdot)$, and an affine map:

$\sigma^{(\ell)}_{W^{(\ell)},b^{(\ell)}}(Z^{(\ell-1)}):=\sigma^{(\ell)}\left(W^{(\ell)}Z^{(\ell-1)}+b^{(\ell)}\right),$
(7)

where $Z^{(\ell-1)}$ is the output from the previous layer, $\ell-1$. The
activation functions, $\sigma^{(\ell)}(\cdot)$, e.g.
$\sigma^{(\ell)}(\cdot)=\max(\cdot,0)$, are critical to non-linear behavior of
the model. Without them, $G_{W,b}$ would be a linear map and, as such, would
be incapable of capturing interaction effects between the inputs. This is true
even if the network has many layers.

##### Deep learning theory

The choice of feed-forward architecture is motivated on theoretical and
empirical grounds. First, feed-forward networks are furnished with a universal
approximation theorem (Hornik et al., 1989). It has further been shown that
deep networks can achieve superior performance versus linear additive models,
such as linear regression, while avoiding the curse of dimensionality (Poggio,
2016).

Martin and Mahoney (2018) show that deep networks are implicitly self-
regularizing and Tishby and Zaslavsky (2015) characterizes the layers as
’statistically decoupling’ the input variables.

There are additionally many recent theoretical developments which characterize
the approximation behavior as a function of network depth, width and sparsity
level (Polson and Rockova, 2018). Recently Bartlett et al. (2017) prove upper
and lower bounds on the expressability of deep feedforward neural network
classifiers with the piecewise linear activation function, such as ReLU
activation functions. They show that the relationship between expressability
and depth is determined by the degree of the activation function.

These bounds are tight for almost the entire range of parameters. Letting $n$
denote the total number of weights, they prove that the VC-dimension is
$\mathcal{O}(nLlog(n))$.

They further showed the effect of network depth on VC-dimension with different
non-linearities: there is no dependence for piecewise-constant, linear
dependence for piecewise-linear, and no more than quadratic dependence for
general piecewise-polynomial. Thus the relationship between expressability and
depth is determined by the degree of the activation function.

There is further ample theoretical evidence to suggest that shallow networks
can’t approximate the class of non-linear functions represented by deep ReLU
networks without blow-up. Telgarsky (2016) shows that there is a ReLU network
with $L$ layers and $U$ units such that any network approximating it with only
$\mathcal{O}(L^{1/3})$ layers must have $\Omega(2^{L^{1/3}})$ units. Mhaskar
et al. (2016) discuss the differences between composition versus additive
models and show that it is possible to approximate higher polynomials much
more efficiently with several hidden layers than a single hidden layer.

###  3.1 PLS Estimation

In the literature, there are multiple types of algorithms for finding the
projections (Manne, 1987). Algorithms differ on whether they estimate the
score matrix $V$ as orthonormal columns or not. The original one proposed by
Wold et al. (1984) uses the conjugate-gradient method (Golub and Van Loan,
2013) to invert matrices. The first PLS projection $p$ and $q$ is found by
maximizing the sample covariance between the $X$ and $Y$ scores

$\displaystyle\operatorname*{maximize\quad}_{p,q}$
$\displaystyle\left(Xp\right)^{\top}\left(Yq\right)$ subject to
$\displaystyle||p||=||q||=1.$

Then the corresponding scores are

$v=Xp,\text{ and }u=Yq.$

We can see from the definition that the directions (loadings) for $Y$ are the
right singular vectors of $X^{\top}Y$ and loadings for $X$ are the left
singular vectors. The next step is to perform regression of $V$ on $U$, namely
$U=V\beta$. The next column of the projection matrix $P$ is found by
calculating the singular vectors of the residual matrices
$(X-vp^{\top})^{\top}(Y-u\beta q^{\top})$.

The final regression problem is solved $Y=UQ=V\beta Q=XP^{\top}\beta Q$. Thus
the PLS estimate is

$\beta_{\mathrm{PLS}}=P^{\top}\beta Q.$

Helland (1990) showed that PLS estimator can be calculated as

$\beta_{\mathrm{PLS}}=R(R^{\top}S_{xx}R)^{-1}R^{\top}S_{xy}$

where $R=(S_{xy},S_{xx}S_{xy},\ldots,S_{xx}^{q-1}S_{xy})$,

$S_{xx}=\dfrac{X^{\top}(I-{\bf 11^{\top}}/n)X}{n-1},$ $S_{xy}=ave(YX),$

where the average is taken over training data, that is
$ave(X)=\frac{1}{n}\sum_{i=1}^{n}\mathbf{x}_{i}$. Once the loadings $(Q,P)$
and scores $(U,V)$ have been estimated, we can then train a deep learner on
$(U,V)$ and then use $U$ to predict $V$ and hence obtain an estimate of $Y$.

To understand why the deep-learner can be decoupled from the estimation of the
projection matrices, we need to establish that the projection directions in
PLS are invariant to the functional form of the regression between score
matrices, up to a constant of proportionality.

We first recall a key asymptotic result from Naik and Tsai (2000) which builds
on an earlier result on OLS regression by Brillinger (2012), namely that PLS
regression is a consistent estimator of the regression coefficient even when
there exists an unknown non-linear relationship between the $\mathbf{y}$ and
$\mathbf{x}$.

######  Assumption 1.

Assume $\mathbf{x}_{1},\dots,\mathbf{x}_{N}$, where
$\mathbf{x}^{\top}=(x_{1},\dots,x_{p})$, are independent normals with mean
$\mu_{x}$ and non-singular covariance matrix $\Sigma_{p\times p}$, and that
the xs are independent of $\bm{\epsilon}_{1},\dots,\bm{\epsilon}_{N}$, where
each $\bm{\epsilon}^{\top}=(\epsilon_{1},\dots,\epsilon_{q})$ has finite
variance $\sigma^{2}$.

######  Assumption 2.

Let $\mathbf{y}=\mathbf{g}+\bm{\epsilon}$, where
$\mathbf{y}^{\top}=(y_{1},\dots,y_{q})$,
$\mathbf{g}^{\top}=(g_{1},\dots,g_{q}),\leavevmode\nobreak\
g_{j}:=g_{j}(z_{j}),z_{j}=\beta_{0j}+\mathbf{x}^{\top}\bm{\beta}_{j},\leavevmode\nobreak\
\beta_{0j}$ is a scalar and $\bm{\beta}_{j}\in\mathbb{R}^{p}$. Assume that
$g(\cdot)$ is a measurable function of $z$ with
$\mathbb{E}[|g_{j}(z_{j})|]<\infty$ and
$\mathbb{E}[|z_{j}g_{j}(z_{j})|]<\infty$ for $j=1,\dots,q$.

######  Theorem 2 (Naik and Tsai (2000)).

Assume the above two assumptions hold. Then the PLS estimator $\hat{B}_{PLS}$
is a strongly consistent estimator of
$B=(\bm{\beta}_{1},\dots,\bm{\beta}_{q})$ up to constants of proportionality.

The theorem is closely related to Stein’s Lemma Stein (2012), which states
that the covariance between the Gaussian regressand and a non-linear Lipschitz
function of the regressors can be written as the covariance between the
regressand and the regressors themselves up to a multiplicative term in the
expectation of the derivative of the non-linear function.

The practical implication of Brillinger’s theorem Brillinger (2012) or Naik &
Tsai theorem is that OLS or PLS respectively generalize to the case
$\mathbb{E}[\mathbf{y}\leavevmode\nobreak\ |\leavevmode\nobreak\ \mathbf{x}]$
is not linear when $\mathbf{x}$ is jointly Gaussian.

For the purpose of modifying the PLS algorithm to use a deep-learner rather
than the linear PLS estimator, we would like to know that the strong
consistency result still holds if we introduce a non-linear regression to
estimate the relation between the scores. Of course, the mere exercise of
reintroducing non-linear regression when Naik & Tsai’s theorem shows that a
linear PLS algorithm will generalize to a non-linear function seems like a
tautology.

However, our immediate motivation, before building a regression model for non-
Gaussian data, is to first show that we needn’t restrict the PLS algorithm to
a linear estimator for the result of Naik & Tsai to hold. Rather, the power of
PLS is in its composability — we can use the x-scores and y-scores estimated
by the PLS algorithm composed with a deep learning estimate of the map between
the scores, rather than needing to substantially overhaul the PLS method to
embed machine learning methods. More precisely, we first seek to generalize
Naik & Tsai theorem by showing that PLS can be composable with another machine
learning method while preserving its consistency property.

In more detail, we show that the strongly consistent estimator,
$\hat{B}_{PLS}$, is invariant under a parameterized non-linear function in the
projected space, $G(\cdot)\equiv G_{\theta}(\cdot)$, s.t.
$u_{j}=G(\beta_{0j}+\mathbf{v}^{\top}\bm{\beta}_{j})$, up to a constant of
proportionality. The x-scores, $\mathbf{v}$, are jointly Gaussian if the
predictors, $\mathbf{x}$, are jointly Gaussian, since the projection,
$P\mathbf{v}$, is a linear transformation.

######  Definition 3.1.

Given the system of equations:

$\displaystyle\mathbf{y}$ $\displaystyle=$ $\displaystyle
Q\mathbf{u}+\bm{\epsilon}_{y},$ $\displaystyle\mathbf{x}$ $\displaystyle=$
$\displaystyle P\mathbf{v}+\bm{\epsilon}_{x},$ $\displaystyle\mathbf{u}$
$\displaystyle=$ $\displaystyle G(\mathbf{v}^{\top}B),$

where $P\in\mathbb{R}^{p\times q}$ and $Q\in\mathbb{R}^{q\times K}$ are
projection matrices referred to as loadings.
$\mathbf{u},\mathbf{v}\in\mathbb{R}^{K}$ are the y-scores and orthonormal
x-scores respectively. $\bm{\epsilon}_{y}\in\mathbb{R}^{q}$ and
$\bm{\epsilon}_{x}\in\mathbb{R}^{p}$ are independent and have finite variance
$\Sigma_{\epsilon_{y}}$ and $\Sigma_{\epsilon_{x}}$ respectively. Let
$w_{j}=\mathbf{v}^{\top}\bm{\beta}_{j}$.

######  Assumption 3.

Assume that $\mathbf{v}$ is Gaussian and is a “good” predictor of $\mathbf{u}$
through some non-linear transformation
$\mathbb{E}[\mathbf{u}\leavevmode\nobreak\ |\leavevmode\nobreak\
\mathbf{v}]=G(\mathbf{v}^{\top}B)$, where $G(\cdot)$ is a measurable function
with $\mathbb{E}[|G_{j}(w_{j})|]<\infty$ and
$\mathbb{E}[|w_{j}G_{j}(w_{j})|]<\infty$ for $j=1,\dots,q$.

######  Theorem 3 (Composability of PLS).

Under this assumption on $G(\cdot)$, embedded in the PLS method, the PLS
estimator, $\hat{B}_{PLS}$, is a consistent estimator of
$B=(\bm{\beta}_{1},\dots,\bm{\beta}_{q})$.

###### Proof.

$\displaystyle\mathbf{y}$ $\displaystyle=$ $\displaystyle
G(\mathbf{v}^{\top}B)Q+\bm{\epsilon},$ $\displaystyle\mathbf{y}$
$\displaystyle=$ $\displaystyle
G((\mathbf{x}-\bm{\epsilon}_{x})P^{\top}\bm{\beta})Q+\bm{\epsilon},$

where $\bm{\epsilon}:=\bm{\epsilon}_{t}^{\top}Q+\bm{\epsilon}_{y}$.

$\displaystyle Cov(y_{j},\mathbf{x})$ $\displaystyle=$ $\displaystyle
Cov(G_{j}(w_{j})Q,\mathbf{x})$ $\displaystyle=$ $\displaystyle
Cov(G_{j}(w_{j}),Q\mathbf{x})$ $\displaystyle=$ $\displaystyle
Cov(w_{j},Q\mathbf{x})\frac{Cov(G_{j}(w_{j}),w_{j})}{Var(w_{j})}$
$\displaystyle=$ $\displaystyle
Cov((\mathbf{x}P^{\top}B)_{j},Q\mathbf{x})\frac{Cov(G_{j}(w_{j}),w_{j})}{Var(w_{j})}$
$\displaystyle=$ $\displaystyle
Cov((\mathbf{x}P^{\top}B)_{j}Q,\mathbf{x})\kappa_{j}$ $\displaystyle=$
$\displaystyle Cov(\mathbf{x},\mathbf{x})(P^{\top}B)_{j}Q\kappa_{j}$

and the PLS estimator
$(\hat{\bm{\beta}}_{PLS})_{j}=\Sigma_{x}^{-1}Cov(y_{j},\mathbf{x})\kappa_{j}$
is a consistent estimator of $P^{\top}\bm{\beta}_{j}Q,\leavevmode\nobreak\
j=1,\dots,q$. ∎

We discuss and speculate on this result.

######  Remark 3.

Since the PLS estimator is invariant under the choice of non-linearity in
functional approximation between the x-scores and y-scores, we can use the
x-scores and projection matrices estimated by the PLS method without needing
to modify them. This provides the justification for composability of PLS with
the deep leaner.

######  Remark 4.

While DPLS is intended to be a general methodology and not limited to Gaussian
predictors, we make no such estimator consistency claim when $\mathbf{x}$ is
not Gaussian. Rather, the stated goal of composability is simply to provide a
better predictor of the Y-scores under the assumption that
$\mathbb{E}[\mathbf{y}\leavevmode\nobreak\ |\leavevmode\nobreak\ \mathbf{x}]$
is non-linear as the linear estimator is no longer “valid” due to the non-
Gaussian predictors.

######  Remark 5.

Under an idealized data generation process, we may use PLS as a baseline for
measuring the performance of DPLS and thus if $\mathbf{x}$ is Gaussian, then
we might expect to characterize any relative gain in performance when using
DPLS in place of PLS under infinite sample sizes. After all, in such a case, a
linear approximation should suffice up to a constant of proportionality.
However, we exercise some caution in over-interpreting the theorem – the
Brillinger and Naik & Tsai result is simply an asymptotic consistency result
and says nothing about whether linear regression is “optimal” when there is a
non-linear relationship between the response and the Gaussian predictors. In
fact the asymptotic covariance matrix used for the least squares estimate of
the regression coefficients depends on the unknown form of the non-linearity
between $\mathbf{y}$ and $\mathbf{x}$. Needless to say, the limitations of an
asymptotic result to sampling error from use of a finite sample size are part
and parcel of any frequentist approach.

######  Remark 6.

One can estimate the weights in the first layer of the network using the PLS
estimator. Put differently, we would expect the regression coefficients in the
PLS method to match the weights in the first layer of the network when
$\mathbf{x}$ is Gaussian, irrespective of the value of the weights and choice
of architecture of the remaining network layers. In practice, the comparison
will not be exact due to sampling noise, but also because of the constant of
proportionality. However, the functional form of $G=G_{W,b}(\cdot)$ is known
and the constant of proportionality can be estimated. When the data is non-
Gaussian, we speculate that the weights in the first can be initialized more
effectively than a random initialization under the assumption that they would
be similar under mild departure from Gaussian assumptions.

Our consistency result permits us to compose the loadings estimated by PLS,
under the original assumption of a linear relationship between the input and
output scores, with a non-linear function $G_{W,b}(\cdot)$ to improve the
predictive performance of the output scores (See Figure 3 for an
illustration). The theorem shows that the composition of a deep-learner with
the projection matrices estimated separately by the PLS algorithm does not
“break” PLS – the asymptotic consistent estimate result of Naik & Tsai still
holds. With this in mind, we propose the following two-step DPLS algorithm
here:

  1. 1.

Assume that the predictors, $X$, are Gaussian and hence by the Naik & Tsai
Theorum, the PLS estimator is a consistent estimator of the coefficients in
$g(\beta_{0}+\beta x_{j})$.

  2. 2.

Relax the assumption of Gaussian predictors: apply a deep-learner to the non-
Gaussian x-scores and y-scores estimated by PLS, estimate a non-linear map to
improve the prediction of the y-scores. Optionally use the PLS coefficients as
initial guesses of the weights in the first layer of $G_{W,b}(\mathbf{v})$.

In this respect, DPLS could be viewed, as an end-to-end deep learning
architecture between $\mathbf{x}$ and $\mathbf{y}$:

$\hat{\mathbf{y}}=Q\circ(\sigma^{(L-1)}_{W^{(L-1)},b^{(L-1)}}\circ\dots\circ\sigma^{(2)}_{W^{(2)},b^{(2)}})(P^{+}\mathbf{x}),$

where $P^{+}=$ in given by the PLS algorithm and is projection of $\mathbf{x}$
on to the first layer, the next layer is non-linear in the scores (with the
weights $(W^{(2)},b^{(2)})$ possibly even estimated by PLS), the penultimate
layer is a linear projection on to space of y-scores (i.e. the activation
function $\sigma^{(L-1)}$ is linear and the final layer is the linear
projection $Q$, estimated by PLS, from the $K$-vector space to the $q-$ vector
space of $\mathbf{y}$.

Figure 3: The y-scores ($K=1$) are shown against the PLS coefficient
multiplied by the orthonormal x-scores to illustrate the non-linearity which
must be captured by the deep-learner when the scores and returns are non-
Gaussian.

More generally, PLS can be used as a layer at any stage of deep learning, not
just the first. For example, in a convolutional neural network,

$Z_{1}=g\left(\sum_{i\in M}\mathbf{x}_{i}\star w_{i}+b_{i}\right),$

where $X\star w+b$ denotes the convolution over the region $M$, with input
$X$, weights $w$ and bias $b$. Then we can add a PLS layer by regressing the
output $Y$ on $Z_{1}$ and perform feature reduction, which results in a CNN-
PLS model.

Our theorem and algorithm is not limited to use with deep learning either. PLS
can be combined with many well-known machine learning models, such as kernel
regression (Rosipal and Trejo (2001), Rosipal et al. (2001)), Gaussian
processes (Gramacy and Lian (2012)), tree models. For example, Higdon et al.
(2008) uses a PCA decomposition of the $Y$-variable before applying a
nonlinear regression method. Of course deep learning is attractive because of
its universal representation theorem and composability with other algebraic
transformations.

###  3.2 PLS Algorithm

One common problem with estimation of linear prediction rules,
$\hat{Y}=X\beta$, on high dimensional datasets is the high variance due to the
typically ill-conditioned nature of the design matrix, $(X^{T}X)$. Shrinkage
methods are used to address this issues and work by biasing coefficient vector
away from directions in which the predictors have low sampling variability— or
equivalently, sway from the ”least important” principal components of $X$.

The Bayesian shrinkage interpretation of Frank and Friedman (1993) provides a
useful perspective for explaining the PLS algorithm. PLS, in fact, relies on a
non-standard shrinkage mechanism which shrinks away from the origin for
certain eigen-directions. The shrinkage factors (a.k.a. scale factors) along
each of the eigen-directions are nonlinear in the response values. They also
depend on the OLS solution.

Polson and Scott Polson and Scott (2010, 2012) provide a general theory of
global-local shrinkage and, in particular, analyze g-prior and horseshoe
shrinkage. Bhadra et al. (2020) build on this work by demonstrating that
horseshoe regularisation is useful for a wide class of machine learning
methods.

We first center and standardise $(Y,X)$. Then, we provide an SVD decomposition
of the sample covariance matrix:

$V=ave(XX^{\top}).$

This finds the eigenvalues $\\{e_{j}^{2}\\}_{j=1}^{p}$ and corresponding
eigenvectors $\\{\bm{v}_{j}\\}_{j=1}^{p}$ arranged in descending order, so we
can write the eigenvector decomposition

$V=\sum_{j=1}^{p}e_{j}^{2}{\bm{v}}_{j}{\bm{v}}_{j}^{\top}.$

This leads to a sequence of regression models
$(\hat{Y}_{0},\dots,\hat{Y}_{R})$ with $\hat{Y}_{0}$ being the response mean
and

$\hat{Y}_{K}=\sum_{k=0}^{K}(ave(w_{k}^{\top}X)/e_{k}^{2})\bm{v}_{k}^{\top}X,\leavevmode\nobreak\
\leavevmode\nobreak\ d\in\\{0,\dots,R\\}.$

where $R$ is the rank of $V$ (number of nonzero $e_{k}^{2}$). Therefore, PLS
finds ”features” $\\{\bm{f}_{k}\\}_{k=0}^{K}$=
$\\{\bm{v}_{k}^{\top}X\\}_{k=0}^{K}$.

The PLS estimator is of the form

$\hat{Y}^{M}=\sum_{j=1}^{R}f_{j}^{PLS}\hat{\alpha}_{j}\bm{v}_{j}^{\top}\bm{x},$

where $f_{j}^{PLS}$ are PLS scale factors for the $K$ top eigenvectors (Frank
and Friedman, 1993) of the form

$\displaystyle f_{j}^{PLS}$
$\displaystyle=\sum_{k=1}^{K}\theta_{k}e_{j}^{2k}\;,\;{\rm
where}\;\theta=w^{-1}\eta,\;\eta_{k}=\sum_{j=1}^{p}\hat{\alpha}_{j}^{2}e_{j}^{2(k+1)}.$

Figure 4: The scale factors are plotted for $j=1,2,3$. If any $f_{j}>1$ then
one can expect supervised learning will lead to different predictions than
unsupervised learning.

As shown in Figure 4, the scale factors provide a diagnostic plot: $f_{j}$. If
any $f_{j}>1$ then one can expect supervised learning (a.k.a. PLS with $Y$’s
influencing the scaling factors) will lead to different predictions than
unsupervised learning. In this sense, PLS is an optimistic procedure in that
the goal is to maximise the explained variability of the output in sample with
the hope of generalising well out-of-sample. For linear estimators, $f_{j}>1$
means that both the bias and the variance are increased. Once the projection
matrices are estimated, a standard stochastic gradient descent algorithm, such
as Adam (Kingma and Ba, 2014), can be used to fit the deep learning model to
predict the y-scores.

###  3.3 Covariate sensitivities

We use a gradient-based technique for ranking the importance of the covariates
in the fitted DPLS model. Note, for avoidance of doubt, that in contrast to
PCA and autoencoders, the DPLS model sensitivities are given in terms of the
observed covariates and not just the latent statistical factors. Our proposed
method is consistent with how coefficients are interpreted in linear
regression, i.e. as _model_ sensitivities and is thus appealing to
practitioners accustomed to use OLS regression. In our approach, the model
sensitivities are the partial derivatives of the fitted model output w.r.t.
each input. This method is consistent with Horel and Giesecke (2019) who
develop statistical tests for the significance of the input variables in a
neural network, using their partial derivatives. The test statistic is based
on a weighted distribution of the square of the neural network partial
derivative, w.r.t. each input. However, the scope of their study is limited to
a single-hidden layer network and treats the asymptotic distribution of the
network as a Gaussian process.

Here we shall consider deep networks with an arbitrary number of layers and do
not rely on asymptotic approximations, which may be limited when the network
has a small number of neurons. Indeed we derive closed form expressions for
the Jacobian and Hessian - the off-diagonals provide sensitivities to
interaction terms. In contrast to Horel and Giesecke (2019), our goal is to
characterize the dispersion of the empirical sensitivity distribution -
confidence intervals of the distributions can be found by Bootstrap sampling
(see, for example, Wood (2005)). Such an approach is only useful in practice
if the variance of these distributions are bounded. Our primary theoretical
result then is to characterize the upper bound on the variance of the partial
derivatives and show that it is bounded for finite weights.

To evaluate fitted model sensitivities analytically in the DPLS method, we
require that the function $\hat{V}=G_{\hat{\theta}}(U)$ is continuous and
differentiable everywhere. Furthermore, for stability of the interpretation,
we shall require that $G_{\hat{\theta}}(\mathbf{x})$ is a Lipschitz
continuous666If Lipschitz continuity is not imposed, then a small change in
one of the input values could result in an undesirable large variation in the
derivative. That is, there is a positive real constant $\kappa$ s.t.
$\forall\mathbf{x}_{1},\mathbf{x}_{2}\in\mathbb{R}^{p}$,
$|G_{\hat{\theta}}(\mathbf{x}_{1})-G_{\hat{\theta}}(\mathbf{x}_{2})|\leq\kappa|\mathbf{x}_{1}-\mathbf{x}_{2}|$.
Such a constraint is necessary for the first derivative to be bounded and
hence amenable to the derivatives, w.r.t. to the inputs, providing
interpretability.

Fortunately, provided that the weights and biases are finite, each semi-affine
function is Lipschitz continuous everywhere. For example, the function
$\tanh(x)$ is continuously differentiable with derivative, $1-\tanh^{2}(x)$,
is globally bounded. With finite weights, the composition of $\tanh(x)$ with
an affine function is also Lipschitz. Clearly ReLU$(x):=\max(\cdot,0)$ is not
continuously differentiable and one must instead use a softplus function, a
smooth approximation to ReLU.

From the definition of the DPLS model, we can write the fitted model in the
compact form

$\hat{Y}=G_{\hat{\theta}}(P^{\top}X)Q.$

Taking first and second derivatives of the response w.r.t. the input variables
gives

$\displaystyle\frac{\partial\hat{Y}}{\partial X}$ $\displaystyle=$
$\displaystyle\frac{\partial G_{\hat{\theta}}}{\partial V}P^{\top}Q$ (8)
$\displaystyle\frac{\partial^{2}\hat{Y}}{\partial X^{2}}$ $\displaystyle=$
$\displaystyle P\frac{\partial^{2}G_{\hat{\theta}}}{\partial V^{2}}P^{\top}Q$
(9)

where $\frac{\partial G_{\hat{\theta}}}{\partial V}$ and
$\frac{\partial^{2}G_{\hat{\theta}}}{\partial V^{2}}$ are the Jacobian and
Hessian matrices of the fitted deep learner w.r.t. to the X-scores, $V$. Thus,
the covariate sensitivities in the DPLS model are linear transformations of
the sensitivities of the fitted deep learner w.r.t. to the latent X-scores.
See Section B, for further details of the functional form of the Jacobian and
Hessian matrix.

###  3.4 DPLS Factor Models

To fix ideas, we begin by presenting a linear PLS factor model and then show
this generalizes to our deep PLS factor model. A PLS factor model is of the
form

$\mathbf{r}^{e}_{t}=\hat{U}Q+\epsilon_{t}=\mathbf{z}_{t-1}P^{+}_{t}B_{t}Q_{t}+\epsilon_{t}=V_{t}(B_{t}Q_{t})+\epsilon_{t}=V_{t}\mathbf{f}_{t}+\epsilon_{t},$
(10)

$P^{+}$ is the Moore-Penrose Inverse of $P$, $B\in\mathbb{R}^{K\times K}$ is a
diagonal matrix with the regression coefficients on the diagonal estimated by
the PLS algorithm.

The latent risk factors, $\mathbf{f}_{t}=B_{t}Q_{t}$, are given by the
projection of the regression coefficients and the factor loadings (the
X-scores), $V_{t}$, are projections of the firm characteristics,
$\mathbf{z}_{t-1}$ onto the K-factors.

In any time period, $t$, we can therefore attribute the expected excess
returns of stock $i$ to the latent factors:

$\bar{r}^{e}_{i,t}:=\mathbb{E}_{t}[r^{e}_{i,t}]:=\mathbb{E}[r^{e}_{i,t}\leavevmode\nobreak\
|\leavevmode\nobreak\
z_{i,t-1}]=\mathbf{v}_{i,t}^{\top}\mathbb{E}_{t}[\mathbf{f}_{t}],\leavevmode\nobreak\
i=1,\dots N_{t}.$ (11)

Similarly, in any period, $t$, we can attribute the conditional variance of
the excess asset returns to the latent factors

$\sigma^{2}_{i,t}:=\mathbb{V}_{t}[r^{e}_{i,t}]:=\mathbb{V}[r^{e}_{i,t}\leavevmode\nobreak\
|\leavevmode\nobreak\
z_{i,t-1}]=\mathbf{v}_{i,t}^{\top}\Sigma_{t}^{f}\mathbf{v}_{i,t}+\tilde{\sigma}^{2}_{i,t},\leavevmode\nobreak\
i=1,\dots N_{t},$ (12)

where $\Sigma^{f}_{t}$ is the covariance matrix of the latent factors at time
$t$ and $\tilde{\sigma}^{2}_{i,t}$ is the variance of the error.

Note that PLS does not explicitly find the explained variance of the returns.
Rather for a target explained covariance between the excess returns and the
characteristics, PLS determines the number of factors $K$.

We now turn to a specific subset of deep fundamental factor models — DPLS
models. DPLS projects the response and predictors below as

$\displaystyle\mathbf{r}^{e}_{t}$
$\displaystyle=U_{t}Q_{t}+\bm{\epsilon}^{y}_{t},$ $\displaystyle Z_{t-1}$
$\displaystyle=V_{t}P_{t}+\bm{\epsilon}^{x}_{t}$

and a feedforward neural network is used to predict the output scores from the
input-scores:

$U_{t}=G_{W_{t},b_{t}}(V_{t})+\mathbf{\epsilon}_{t}.$ (13)

The primary benefit of using deep learning, is that it is able to
automatically capture non-linear effects and, in particular, is robust to
outliers.

The DPLS model gives the risk premia $V_{t}$ on the latent risk factors by a
Taylor expansion of the deep learner, $g:=G_{W,b}(V_{t})$, about the origin
$\mathbf{0}$:

$\displaystyle\mathbf{r}^{e}_{t}$ $\displaystyle=$
$\displaystyle\hat{U}_{t}Q_{t}+\epsilon_{t}$ (14) $\displaystyle=$
$\displaystyle g(V_{t})Q_{t}+\epsilon_{t}=g(\mathbf{0})Q_{t}+V_{t}\nabla
g(\mathbf{0})Q_{t}+\frac{1}{2}V_{t}\nabla^{2}g(\mathbf{0})V_{t}Q_{t}$ (15)
$\displaystyle=$
$\displaystyle\alpha_{t}+V_{t}\mathbf{f}_{t}+\frac{1}{2}V_{t}H_{0}V_{t}Q_{t},$
(16)

where the latent risk factors are given by $\mathbf{f}_{t}=\nabla
g(\mathbf{0})Q_{t}$ and the quadratic term, defining the non-linear factor
structure, requires the Hessian, $H_{0}$, evaluated at the origin. This
Hessian contains the interaction effects w.r.t. to the risk premia (x-scores).

####  3.4.1 State-space DPLS fundamental Factor Models

For completeness, we can explore other variants of the DPLS fundamental factor
model above – specifically it can be written in a more general form which
includes higher order lag terms up to maximum lag $r$. Specifically we posit
the non-linear state-space model:

$\displaystyle\mathbf{r}^{e}_{t}$ $\displaystyle=$ $\displaystyle
U_{t}Q_{t}+\bm{\epsilon}^{y}_{t},$ $\displaystyle Z_{t-1}$ $\displaystyle=$
$\displaystyle V_{t}P_{t}+\bm{\epsilon}^{x}_{t}$

with the score state-vectors defined as

$U_{t}=\begin{bmatrix}U_{t}\\\ \cdots\\\ U_{t-r}\end{bmatrix},\qquad
V_{t}=\begin{bmatrix}V_{t}\\\ \cdots\\\ V_{t-r}\end{bmatrix},$

where $Q_{t}\in\mathbb{R}^{q\times Kr}$ and $P_{t}\in\mathbb{R}^{p\times Kr}$.
Note, as a special case, that the vector autoregressive form (a.k.a. semi-
supervised learning) of the DPLS model (i.e. $X=Y$) is the system of equations

$\displaystyle X_{t}$ $\displaystyle=$ $\displaystyle
U_{t}Q_{t}+\bm{\epsilon}^{x}_{t},$ $\displaystyle U_{t}$ $\displaystyle=$
$\displaystyle G_{\bm{\theta}_{t}}(U_{t-1})+\mathbf{e}_{t}$

where $\mathbf{e}_{t}^{\top}=(\mathbf{\epsilon}_{t},0,\dots,0)$.

If the projection matrix and parameters, $\bm{\theta}$ are further fixed in
time then we see that the DPLS factor model is in fact a non-linear Dynamic
Factor Model (DFM) (Stock and Watson, 2010) with projections on to scores:

$\displaystyle X_{t}$ $\displaystyle=$ $\displaystyle
U_{t}Q+\bm{\epsilon}^{x}_{t},$ $\displaystyle U_{t}$ $\displaystyle=$
$\displaystyle G_{\bm{\theta}}(U_{t-1})+\mathbf{e}_{t},$

where $U_{t}$ is the latent state variable of factors (a.k.a. score state
vector). In the Stock and Watson (2010) model, $Q$ is a linear projection onto
lagged latent factors and $G_{\bm{\theta}}(\cdot)$ is replaced with a static
state-space update matrix. Their model requires stationarity. One important
advantage of our DPLS model, aside from the ability to capture non-linearity
in the latent state update, is that $G_{\bm{\theta}}$ can be an advanced
recurrent neural network architecture such as LSTMs or GRUs, thus relaxing the
requirement for stationarity (Dixon, 2022).

##  4 Experiments

This section presents the application of our DPLS factor model to all stocks
in the Russell 1000 index, which represents the top 1000 companies by market
capitalization in the United States. The historical data covers an
approximately 30 year period of monthly updates between January 1989 to
November 2018, with a coverage universe of 3290 stocks. In every month, each
stock includes 49 BARRA model factors — 19 of which are fundamental factors
and the remainder are GICS sector dummy variables (see Section A for a
description of the factors).

Note that the BARRA factor model includes many more explanatory variables than
used in our experiments below, but the purpose, here, is to illustrate the
application of our framework to a larger set of factors. Equally, the BARRA
factors are designed to explain variance in asset returns rather than capture
alpha, however, our results focus both on expected asset returns and explained
variance of asset returns based on conditional exposures to latent risk
factors. Note, for avoidance of doubt, that our goal in this section is not to
shed further light on the “factor zoo problem”, but simply to advocate for the
use of DPLS as a methodology for factor modeling based on its ability to
capture outliers and attribute prediction error to conditional exposures to
latent risk factors.

Our experiments broadly focus on three aspects of the DPLS factor model:

  1. 1.

Ability to capture outliers in monthly excess returns;

  2. 2.

Ability to generate higher performing managed portfolios under a risk adjusted
return metric; and

  3. 3.

Ability to explain asset returns and risk via risk premia.

##### Data cleaning:

All stocks with missing factor exposures are removed from the estimation
universe for the date of the missing factors only. To avoid excessive turn-
over in the estimation universe over each consecutive period, we include all
dropped symbols in the index over the 12 consecutive monthly periods. Further
details of the data preparation, including sanitization to avoid violating
licensing agreements, are provided in the documentation for the source code
repository referenced on the first page.

##### Rescaling:

Training and testing is alternated in each period. For example, in the first
historical month of the data, the model is fitted to the factor exposures and
monthly excess monthly returns over the next period. To avoid look-ahead bias,
all factors in each training period are normalized using only the training
data, i.e. the resulting mean and standard deviation of each factor across all
stocks is 0 and 1 respectively. The same normalization parameters are then
used to rescale the test set. This avoids introducing a bias into the test set
while also avoiding incorporating knowledge of the test data when rescaling
the training data. Of course, the test set is not perfectly normalized.

##### Hyper-parameter tuning:

Cross-validation is performed using this cross-sectional training data, with
approximately 1000 symbols. Once the model is fitted and tuned, we then apply
the model to the factor exposures in the next period, t+1, to predict the
excess monthly returns over [t+1, t+2]. The process is repeated over all
periods. Note, for computational reasons, we can avoid cross-validation over
every period and instead stride the cross-validation, every other say 10
periods, relying on the optimal hyper-parameters from the last cross-
validation periods for all subsequent intermediate periods. In practice we
find that performing cross-validation in the first period only is adequate
with the exception of the number of latent risk factors, $K$.

Three fold cross-validation over $\\{50,100,200\\}$ hidden units per layer,
$\lambda_{1}\in\\{0,10^{-3},10^{-2},10^{-1}\\}$ and $\\{1,2,3\\}$ hidden
layers is performed in the same time period. We find that the optimal
architecture has two hidden layers, 100 units per layer, no $L_{1}$
regularization, and softplus (i.e. smooth ReLU) activation. For an
interpretable model, $C^{2}(\mathbb{R})$ continuous activation functions are
required to provide sufficient smoothness (with 50 hidden units per layer
being optimal).

All experiments are performed using `Tensorflow` (Abadi et al., 2016) to
implement the deep feed-forward network and the `pls` R package (Liland et
al., 2021) for PLS. The LASSO model is implemented using the `glmnet` R
package (Simon et al., 2011). The amount of $L_{1}$ regularization in the
LASSO model is tuned every period777The LASSO model uses cross-validation to
optimize the $L_{1}$ regularization parameter and iterative fitting, with 50
alphas, along the regularization path.. The PCA results are generated using
the asymptotic PCA method implemented in the `covFactorModel` R
package888https://rdrr.io/github/dppalomar/covFactorModel.

###  4.1 Russell 1000 stocks

We begin by analyzing the predicted excess assets returns of the models over
all constituents of the Russell 1000 index in each period999Note, for
avoidance of doubt, that we do not construct index tracking portfolios – all
stocks in the index are equally weighted in our analysis., with the goal of
identifying the ability of DPLS to capture outliers in the excess stock
returns. Figure 5 compares the in-sample and out-of-sample performance of
DPLS, PLS, deep feedforward network (NN), and LASSO regression using the
$L_{\infty}$ norm of stock return errors over the coverage universe of stocks
listed in the Russell 1000 index at each monthly period. This norm measures
the largest absolute error across all assets, i.e. the error associated with
the most mis-predicted excess asset return in each period. The time averaged
$L_{\infty}$ norms, over all periods, is shown in parentheses.

We observe the ability of DPLS to best capture outliers, with the $L_{\infty}$
norm of the error in the NN and DPLS regression respectively being on average
between 8% and 12% smaller than the other methods out-of-sample. We also
observe that the amount of variance between the in-sample and out-of-sample
performance is most notable using NN regression. This is consistency with the
high degree of parameters exhibited by the full neural network architecture
relative to the training data size. This, in fact, highlights a substantial
defect of using deep learning for period-by-period cross-sectional regression
— the amount of training data required for deep learning required isn’t
compatible with factor modeling cross-sectional datasets. On the other hand,
if one simply trains across all stocks and all periods, then the model loses
its dynamic nature.

We also observe periods in the out-of-sample test set where PLS exhibits the
largest errors and DPLS appears to capture these. The combination of a lower
time averaged $L_{\infty}$ error and fewer apparent outliers not captured by
the PLS, provides evidence of the superiority of DPLS over PLS in capturing
outliers.

(a) $L_{\infty}$ in-sample error

(b) $L_{\infty}$ out-of-sample error

Figure 5: The (a) in-sample and (b) out-of-sample stock returns error, under
the $L_{\infty}$ norm, is compared between DPLS, PLS, NN, and LASSO regression
applied to a coverage universe of 3290 stocks from the Russell 1000 index over
the period from January 1990 to November 2018. The time averaged $L_{\infty}$
error is shown in parentheses.

Figure 6 compares the in-sample and out-of-sample performance of DPLS, NN, and
LASSO regression using the $R^{2}$ of monthly asset returns. The time averaged
$R^{2}$ over all periods is shown in parentheses. Overall, we observe that
DPLS regression exhibits the best out-of-sample performance. We also find
further evidence that the NN is over-fitting, on account of the low out-of-
sample $R^{2}$.

(a) In sample $R^{2}$

(b) Out-of-sample $R^{2}$

Figure 6: The (a) in-sample and (b) out-of-sample stock returns error, under
the $R^{2}$ metric, is compared between DPLS, PLS, NN, and LASSO regression
applied to a coverage universe of 3290 stocks from the Russell 1000 index over
the period from January 1990 to November 2018. The time averaged $R^{2}$ error
is shown in parentheses.

Figure 7 compares the time averaged factor model sensitivities, evaluated at
$Z_{t-1}=\mathbf{0}$, over the thirty year period using (a) full NN and (b)
DPLS regression. For each method, the sensitivities are each sorted in
descending order from top to bottom by their mean value over all periods. See
Table 3 in Appendix A for a description of the factors.

|
---|---
(a) NN | (b) DPLS
Figure 7: The time averaged factor model sensitivities over the thirty year
period using (a) NN and (b) DPLS regression applied to the constituents of the
Russell 1000 index in each time period. Market Cap (MC), trading activity
(TrA), and Earnings Volatility to Total Assets (EaV/TA) are ranked highly as
positive factors by each model. Additionally, both models agree on Total
Assets (TA), Book-to-Price (B/P), and Cash flow Volatility to Total Assets
(CFV/TA) as the highest ranked negative factors contributing to asset returns.
See Table 3 for a description of the factors.

All three methods share important consistencies. For example, Market Cap (MC),
Trading Activity (TrA), and Earnings Volatility to Total Assets (EaV/TA) rank
highly as positive factors. Conversely, Total Assets (TA), Book-to-Price
(B/P), and Cash flow Volatility to Total Assets (CFV/TA) rank highly as
negative factors. Both models also agree that expected returns are, on
average, largely independent of Cash Flow to Price (CF/P).

However, there are also significant discrepancies between the methods. For
example EBIDTA to Earnings Volatility (EB/EV) ranks highly as a positive
factor in the NN model model, yet is weakly positive in the DPLS model. As
another example, Total Assets Growth (TAG) is ranked as one of the most
negative factors in the DPLS model, yet the NN model places it as weakly
negative. The NN model is observed to exhibit sensitivities which are up to
approximately 50% larger than in the DPLS model, suggesting that the higher
degree of parameterization leads to factor over-sensitivity. The NN
sensitivities are also, on average, skewed negative. Another apparent
difference is the effect of dimension reduction in the DPLS model – there are
many more non-zero fundamental factors in the NN model, whereas the mid-
section of the DPLS plot is comparatively thin, indicating fewer material
factors on average.

(a) Temporal distribution of NN factor sensitivities

(a) Temporal distribution of DPLS factor sensitivities

Figure 8: The temporal distribution of factor model sensitivities over the
thirty year period using (a) NN and (b) DPLS regression applied to the
constituents of the Russell 1000 index in each time period. The medians are
shown by the horizontal solid lines. We observe that Cash flow to Price and
sub-industry group 9 carry the most uncertainty - the inter-quartile range is
the largest in both models. DPLS also collapses a number of the sector and
sub-industry group variables closely around zero.

The time averages provide a limited view on the model sensitivities. We
additionally plot the temporal distribution of the factor sensitivites. We
observe that Cash flow to Price and sub-industry group 9 carry the most
uncertainty - the inter-quartile range is the largest in both models. DPLS
also collapses a number of the sector and sub-industry group variables closely
around zero.

We can also study the extent to which the most dominant factors persist over
multiple periods. Figure 9 compares the sensitivity of the predicted excess
asset returns to the fundamental factors over the two year period between
October 2007 and October 2009, covering the Lehman Brothers financial crash.
We find consistency between the factor sensitivities and the September/October
2008 Lehman collapse period. In particular, DPLS is able to identify a number
of prominent factors with large positive sensitivities in this period. In
contrast, the NN model only identify one strongly positive sensitivity. We
also see more variability from period to period in the NN model, whereas the
DPLS exhibits more observations where the factor dominance persists over
multiple periods.

|
---|---
(a) NN | (b) DPLS
Figure 9: The dynamical model sensitivities over a two year period between
October 2007 and October 2009 covering the Lehman Brothers financial crash.
The factors are ranked by their time averaged sensitivities from top to bottom
in descending order.

Further insight can be gained by analyzing the interaction effects and
quadratic terms in the NN and DPLS models, as shown in Figure 10. These
interaction terms are time averaged over the thirty year period using (a) NN
and (b) DPLS regression applied to all stocks in the Russell 1000 index. $XX-
YY$ denote pairwise interaction effects between factors XX and factors YY.

In contrast to the sensitivities, the most important interaction terms are
different across each model. A common pattern is for a model to find a
fundamental ratio across multiple sectors and industry groups. For example,
the DPLS model finds the quadratic Book-to-Price ratio positively dominant,
but also finds Book-to-Price across SubIndustry 5 in addition to Industry
Groups 0 and 3. The DPLS model also finds Cashflow/Price strongly negatively
dominant across Industry 4, SubIndustry 7, and Sectors 7 and 9. The NN model
finds Earnings Volatility/Total Assets (EaV/TA), Cash Flow Volatility/Total
Assets (CFV/TA), Total Asset Growth (TAG), Total assets (TA), and Market Cap
(MC) positively dominant in specific SubIndustries, Sectors, and Industry
Groups. Whereas Dividends (DIV) are exclusively associated with negative
effects in specific SubIndustries, Sectors, and Industry Groups. We also note
that DPLS isolates a few particularly dominant interaction pairs - the
interaction of Book-to-Price and Earnings-to-Price is a strong negative effect
and the Industry 4-Sector 0 pair is strongly positive.

The top interaction effects in the NN model are also an order of magnitude
larger than in the DPLS model, suggesting that the higher degree of
parameterization renders the NN interaction effects as over-sensitive.

|
---|---
(a) NN | (b) DPLS
Figure 10: The most important pairwise interaction effects, time averaged over
a thirty year period using all stocks.

The temporal distribution of interaction effects can shed light on the degree
of variability in interaction effects over time. Figure 11 compares the NN and
DPLS interaction terms, listed in alphabetic order from left to right. In
contrast to the NN model, most of the DPLS interaction effects are tightly
distributed with a few notable exceptions - Book-to-Price$-$Earnings/Price
(B/P-E/V), Cash flow to price$-$Industry 4 (CF/P-I4) and cash flow to
price$-$SubIndustry 5 (CF/P-SI5).

The NN interaction effects, on the other hand, are much more widely
distributed and it more difficult, if not impossible, to draw meaningful
conclusions from the time averaged ranking of the interaction effects, given
there is so much variability.

(a) Temporal distribution of NN factor interactions

(a) Temporal distribution of DPLS factor interactions

Figure 11: The temporal distribution of factor model interactions over the
thirty year period using (a) NN and (b) DPLS regression applied to the
constituents of the Russell 1000 index in each time period. The medians are
shown by the horizontal solid lines.

Here we observe much less consistency between the NN and the DPLS model. The
top interaction effects in the NN model are also an order of magnitude larger
than in the DPLS model, suggesting that the higher degree of parameterization
renders the NN interaction effects as over-sensitive.

|
---|---
(a) NN | (b) DPLS
Figure 12: The interaction effects (a) of the NN model over a four year period
between October 2007 and October 2011 and (b) of DPLS over a two year period
between October 2007 and 2009, both periods covering the Lehman Brothers
financial crash. The factors are ranked by their time averaged interactions
from top to bottom in descending order.

Figure 12 compares the most important interaction and quadratic effects over a
multi-year period covering the Lehman Brothers financial crash. The DPLS
interaction terms in October of 2008 dominant the other dates and largely
govern the overall time averaged ranking. In contrast, the NN exhibits
multiple outlier dates and more severe swings between an interaction effect
being positive dominant and negative dominant from month to month. For
example, in Figure 12(a), the interaction effects in the first and third row
swing between dark red and dark blue in a relatively short period of time.

###  4.2 Portfolio analysis

We now evaluate the utility of the DPLS factor model as a predictive signal
for constructing portfolios.

To evaluate the relative portfolio performance under a DPLS factor model
versus a PLS, LASSO or NN factor model, we construct an equally weighted,
long-only, portfolio of $n$ stocks with the highest predicted monthly returns
in each period. The portfolio is reconstructed each period, always remaining
equal but varying in composition. This is repeated over the 30 year period in
the data. We then estimate the information ratios from the mean and volatility
of the excess monthly portfolio returns, using the Russell 1000 index as the
benchmark.

This exercise is merely to assess the relative merits of DPLS factor models
versus other methods as is not intended as an exhaustive study of the
implications of using DPLS factor models fro portfolio construction. One
important caveat is that the information ratio does not account for trading
costs, which could erode performance due to high turnover. Additionally we do
not consider more complex constraints such as sector balancing of the managed
portfolios.

Figure 13: The information ratios of a portfolio selection strategy which
selects the $n$ stocks from the universe with the highest predicted monthly
returns. The information ratios are evaluated for various equally weighted,
long-only, portfolios whose size is shown by the x-axis.

Figure 13 compares the information ratios using DPLS, NN, and LASSO regression
to identify the stocks with the highest predicted returns. The information
ratios are evaluated for equally weighted portfolios of various sizes. Also
shown, for control, are randomly selected portfolios, without the use of a
predictive signal. The mean information ratio for each model, across all
portfolios, is shown in parentheses. We observe that the information ratio of
the portfolio returns, using DPLS regression, is approximately 1.2x greater
than NN regression. NN regression, despite exhibiting substantially lower out-
of-sample error than LASSO regression results in an inferior information ratio
than LASSO, suggesting that stock selection based on more accurate predicted
returns does not directly translate to higher information ratios on account of
the stock return covariance.

We also observe that the information ratio of the baseline random portfolio is
small, but not negligible, suggesting sampling bias and estimation universe
modification have a small effect.

Figure 14 shows the time averaged factor tilts of equally weighted portfolios
constructed from the predicted top performing 50 stocks, in each monthly
period, over all time periods.

Most of the factor tilts are consistent across the models, however, some are
not. For example, Market Cap (MC), Total Asset Growth (TAG), and Dividends
(Div). LASSO and NN regression overweight the TAG. LASSO also overweights
Market Cap. DPLS (95%) and NN regression overweight Dividends. For
completeness, the changes in factor tilts over each period are provided in
Section D.1.

Figure 14: The rescaled factors tilts of equally weighted portfolios
constructed from the predicted top performing 50 stocks and then averaged over
time.

Figure 18 shows the rescaled sector tilts of equally weighted portfolios
constructed from the predicted top performing 50 stocks, and then averaged
over all periods. The sectors are ranked by their time averaged ratios, but
their tilts vary each month due to portfolios turn-over. WN denotes a random
equally weighted portfolio.

Figure 15: The rescaled sector tilts of equally weighted portfolios
constructed from the predicted top performing 50 stocks and then averaged over
time.

Utilities and Information Technology are the most dominant sectors, i.e. with
the highest time averaged representation. In contrast, Consumer Discretionary
and Health Care are the least representative sectors shown. Note that the
least representative sector, Real Estate, has been excluded. The sector tilts
across DPLS, NN, and LASSO regression are found to generally concur. Notable
differences include DPLS (95%) favoring portfolios with a higher Information
Technology and Consumer Staples tilt and lower Industrial tilt. For
completeness, the changes in sector tilts over each period are provided in
Section D.2.

###  4.3 Risk factor exposures

We estimate the K-factor DPLS model for various $K$ using the in-sample total
$R^{2}$:

$R^{2,in}_{total}=1-\frac{\sum^{N_{t},T}_{i=1,t=1}(r_{i,t}-\hat{r}_{i,t}(Z_{t-1};\hat{\theta}_{t}))^{2}}{\sum^{N_{t},T}_{i=1,t=1}r^{2}_{i,t}},$

and the out-of-sample total $R^{2}$:

$R^{2,out}_{total}=1-\frac{\sum^{N_{t},T-1}_{i=1,t=1}(r_{i,t+1}-\hat{r}_{i,t+1}(Z_{t};\hat{\theta}_{t}))^{2}}{\sum^{N_{t},T-1}_{i=1,t=1}r^{2}_{i,t+1}},$

over all periods. These total $R^{2}s$ measure the extent to which the latent
risk factors capture the realized riskiness in the panel of individual stocks
either in-sample or out-of-sample respectively. It can be interpreted as the
amount of explained variance of the asset returns. If the intercept is set to
zero in the DPLS model, each of these $R^{2}s$ represent the model’s ability
to describe risk compensation solely through exposure to the latent risk
factors. Otherwise, they represent a risk-free compensation and, for each $K$,
incremental risk compensation adjustments.

Table 1 compares the performance of three different statistical factor models,
both in-sample and out-of-sample, across all $N_{t}$ stocks in every period
$t$. The motivation for showing the in-sample and out-of-sample performance is
to assess the bias-variance trade-off. In particular, a drop in performance
from in-sample to out-of-sample is indicative of over-fitting whereas poor
performance in-sample suggests under fitting. The table also compares the
performance of an equally weighted portfolio which is long the top ten
performing stocks, as predicted by each model. In that case the sum of $i$
collapses as there is just one sum of $t$ for the portfolio returns.

The first model we consider is the asymptotic PCA model (APCA) which is fitted
to panel data rather than period-by-period. This model is solely based on
excess returns and ignores the factor data. Hence, the model can not be used
for prediction. When applied to the cross-section of stocks, we observe a
monotonic increase in the in-sample total $R^{2}$ with increasing $K$. The
performance based on the total $R^{2}$ of the portfolio returns is very high
and increasing factors brings the in-sample performance to an almost perfect
level. However we note that only a few statistical factors in the PCA already
substantially captures the portfolio risk.

The PLS and DPLS are different to APCA in two important ways: (i) neither PLS
nor DPLS adjust for heteroschedasticity as does APCA; and (ii) APCA fits to
panel data whereas PLS ad DPLS fit to period-by-period cross-sectional data.
Hence in PCA, the latent factors are the same across all periods, whereas in
PLS and DPLS, the number of factors vary in each period.

The methodology for selecting the number of latent risk factors is also more
complex in PLS and DPLS as both are supervised learning methods. In each
period, $t$, the optimal number of factors, $K^{*}_{t}$, is chosen for each
method using all stock data according to the hyper-parameter tuning details
described earlier in this section. Each model is then trained with the optimal
number of factors in each period. This is important, for two reasons: (i) the
comparison between PLS and DPLS must show the attribution of the model error
to factors given that the model fits the data well otherwise $K^{*}$ conflates
quality of model fit with prediction attribution against the number of latent
factors given the best fit; and (ii) it becomes computationally intractable to
retrain the DPLS model in each every period over all values of $K$.

Once $K^{*}_{t}$ is estimated, we predict the asset returns using
$K_{t}\in\\{1,\dots,K^{*}_{t}\\}$ — both PLS and DPLS quite naturally perform
this since averaging over the first $K_{t}$ columns of $\hat{U}_{t}$ is a step
in the prediction. However, it is not possible to predict with more than
$K^{*}_{t}$ factors. This procedure is repeated for all periods and the
$R^{2}_{total}(K)$ is shown in the table for in-sample and out-of-sample
forecasting. Note that, because $K^{*}_{t}$ varies between each period, larger
values of $K$ will exhibit fewer entries in the $R^{2}_{total}$ formula prior
to summing over time. For example, suppose there are two periods only, and
$K^{*}_{1}=1$ and $K^{*}_{2}=2$, then there will be two error estimations at
$K=1$ but only one at $K=2$ in the sum. We choose $K_{max}=10$ heuristically
in the table based on a minimum sample size of 100 over all stocks and periods
at $K=K_{max}$ in order to avoid statistical noise dominating the estimations.
In other words, there must be at least 100 error terms to be “double summed”
over $i$ and $t$ in the total $R^{2}$.

Across all stocks, the PLS and DPLS model’s in-sample performance is
substantially higher than PCA and monotonically increases with K, capturing as
much as 25% of the variance of the asset’s excess returns. At first, this
result seems counter intuitive given that PCA optimizes explained variance
whereas PLS optimizes explained covariance between $X$ and $Y$. However, the
tuning of $K_{t}$ is by out-of-sample MSE. Moreover, the neural network in
DPLS is fitted by minimizing a MSE based loss function. Both PLS and DPLS
perform comparably well in-sample across all stocks.

The out-of-sample performance of PLS starts to sharply deteriorate after 4
factors whereas DPLS peaks at 10 factors with as much as 15% explained
variance and almost monotonically grows with the number of factors. Note that
monotonicity of MSE in the number of factors is only guaranteed in any one
period — over multiple periods, there is no reason why the monotonicity should
be preserved. It just so happens, here, that $K^{*}_{t}$ tends to be close to
10.

Applied to the portfolio, PLS and DPLS perform poorly in-sample compared to
PCA, however, DPLS is shown to outperform PLS significantly out-of-sample.
Out-of-sample explained variances of portfolio returns are only available for
$K=1$ using PLS due to too few samples for the total $R^{2}$ calculation. In
other words, there are too few periods where $K^{*}_{t}>1$ to generate
sufficient terms to evaluate the $R^{2}$ over.

A key observation is that the variance of asset returns is typically explained
by a few latent risk factors in the PLS model whereas the DPLS model
attributes the variance to a greater number of factors. One plausible
explanation for this difference might be model reduction in a linear model—
for any stock in a given period, a linear model is characterized by its
intercept and slope. Thus a small number of systemic factors are needed to
capture the cross-sectional variability of the intercept and slope. In
contrast, the DPLS model is not functionally equivalent nor can it be reduced
to such a low dimensional set of risk factors due to the non-linearity. Thus
one might arrive at different conclusions about the number of important
systemic risk factors driving returns based on the inclusion of non-linearity
or not.

|  $K$=1 |  $K$=2 |  $K$=3 |  $K$=4
---|---|---|---|---
| PCA | PLS | DPLS | PCA | PLS | DPLS | PCA | PLS | DPLS | PCA | PLS | DPLS
Stocks (in-sample) | 0.0370 | 0.2182 | 0.1951 | 0.0589 | 0.2702 | 0.2280 | 0.0788 | 0.2887 | 0.2355 | 0.0976 | 0.2979 | 0.2382
Stocks (out-of-sample) | NA | 0.0858 | 0.0571 | NA | 0.0873 | 0.0735 | NA | 0.0790 | 0.0978 | NA | 0.0838 | 0.0971
Portfolio (in-sample) | 0.6759 | 0.8993 | 0.7434 | 0.8301 | 0.7710 | 0.7630 | 0.8961 | 0.8303 | 0.7579 | 0.9132 | 0.8038 | 0.7351
Portfolio (out-of-sample) | NA | 0.5070 | 0.5272 | NA | NA | 0.6097 | NA | NA | 0.5071 | NA | NA | 0.4656
|  $K$=5 |  $K$=6 |  $K$=7 |  $K$=8
| PCA | PLS | DPLS | PCA | PLS | DPLS | PCA | PLS | DPLS | PCA | PLS | DPLS
Stocks (in-sample) | 0.1132 | 0.3043 | 0.2320 | 0.1757 | 0.3108 | 0.2355 | 0.2697 | 0.3114 | 0.2453 | 0.3470 | 0.3114 | 0.2523
Stocks (out-of-sample) | NA | 0.0489 | 0.1082 | NA | 0.0401 | 0.1240 | NA | 0.0077 | 0.1395 | NA | -0.0039 | 0.1441
Portfolio (in-sample) | 0.9359 | 0.8695 | 0.7357 | 0.9718 | 0.8878 | 0.7666 | 0.9653 | 0.8385 | 0.7757 | 0.9691 | 0.8301 | 0.7889
Portfolio (out-of-sample) | NA | NA | 0.5496 | NA | NA | 0.5807 | NA | NA | 0.6281 | NA | NA | 0.6517

Table 1: A comparison of the in-sample and out-of-sample total $R^{2}$ for
individual stocks and managed equally weighted long-only portfolios for
various numbers of latent risk factors across different statistical factor
models. DPLS is shown to yield higher out-of-sample performance compared to
PLS.

In addition to understanding how the latent risk factors explain the asset
returns and portfolio risk, we can use DPLS to explain the effect of non-
linearity on the in-sample model prediction. After all, the non-linear factor
structure is learned using the neural network composed with the PLS loadings.

Figure 16 shows the excess monthly returns of an equally weighted, long-only,
portfolio of top performing assets, as predicted by the DPLS factor model. The
figures shows the attribution of the in-sample predicted excess monthly
portfolio returns by the top three latent risk factors, quadratic term, higher
order terms (H.O.T.), intercept and remaining risk factors. In any period, the
quadratic term is found by a Taylor expansion about zero in the DPLS model.
The higher order terms are found by subtracting the first two terms of the
Taylor expansion from the DPLS prediction. The numbers in the parentheses
denote the time averaged component values.

In periods of market dislocation, as exemplified by the collapse of Lehman
Brothers in October 2008, we observe substantial quadratic effects (red) in
the construction of the signal which highlight their importance in capturing
the outliers. In other periods, where the volatility is lower, the effect of
non-linearity is much more marginal.

Figure 16: The attribution of the in-sample predicted excess monthly portfolio
returns using the DPLS factor model. The predicted returns are attributed to
the top three latent risk factors, quadratic term, higher order terms
(H.O.T.), intercept and remaining risk factors. The numbers in the parentheses
denote the time averaged component values.

###  4.4 Computational considerations

For completeness, we also comment on the computational advantages of using
DPLS compared to NNs. Of course, PLS and NN require substantial computation to
train in contrast with PLS or LASSO.

Table 2 below records the prediction performance and training times of the
various methods, each cross-validated, for one period in the factor data. We
see that the DPLS training times vary from approximately 2x up to 20x faster
compared to a deep neural network and the number of parameters exceeds in one
case more than 10x in reduction, providing substantial memory savings.

Model Type | #Predictors | #Neurons per Layer | #Parameters | MSE (in-sample) | MSE (out-of-sample) |  $L_{\inf}$ Error (in-sample) |  $L_{\inf}$ Error (out-of-sample) | Training Time (s)
---|---|---|---|---|---|---|---|---
DPLS (50%) | 14 | 100 | 11701 | 0.022 | 0.033 | 0.908 | 0.963 | 3183.806
DPLS (80%) | 28 | 50 | 4051 | 0.021 | 0.030 | 0.909 | 0.971 | 833.72
DPLS (95%) | 37 | 100 | 14001 | 0.018 | 0.037 | 0.798 | 0.978 | 9870.077
NN | 49 | 200 | 50401 | 0.019 | 0.028 | 0.798 | 1.065 | 19452.267

Table 2: The prediction performance and training times of the various methods,
each cross-validated, for one period in the factor data.

##  5 Summary

This paper introduces a new approach to conditional latent factor modeling
which combine PLS and deep learning to predict and explain stock returns. DPLS
projects firm characteristics on to a smaller set of systemic dynamic risk
factors by maximizing the covariance between characteristics and asset
returns. Unlike other deep learning based factor models, DPLS factor models
corresponds to a non-linear stochastic discount factor. Additionally the about
of no-arbitrage can be estimated in the model over multiple periods and, if
needed, the intercepts can be set to zero to enforce static no-arbitrage. DPLS
is shown to identify the non-linear latent factor structure — the convexity
effects in the latent factor exposures can be quantified and used to estimate
the extent to which linear risk premia fall short of explaining the risk
reward trade-off. Further more DPLS is shown to exhibit superior out-of-sample
performance than OLS, LASSO, NNs, and PLS.

While DPLS is in its infancy, we comment on how it might gain adoption by
portfolio managers. In our view, model risk is mitigated by running many
models and establishing a consensus - this is a long established approach used
in metereological forecasting for example.

The best use-case for the DPLS factor model is likely running it side-by-side
a linear factor model and asking how the weights in the portfolios change
under the inclusion of non-linearity. The non-linearity is likely to be a
source of nuisance or a welcomed directional adjustment by portfolio managers
as their exposure limits are adjusted up or down commensurately.

## References

  * Abadi et al. (2016) Abadi, M., P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, et al. (2016).  TensorFlow: A System for Large-scale Machine Learning.  In Proceedings of the 12th USENIX Conference on Operating Systems Design and Implementation, OSDI’16, pp. 265–283.
  * Algaba et al. (2021) Algaba, A., S. Borms, K. Boudt, and B. Verbeken (2021).  Daily news sentiment and monthly surveys: A mixed-frequency dynamic factor model for nowcasting consumer confidence.  International Journal of Forecasting.
  * Bai and Ng (2002) Bai, J. and S. Ng (2002).  Determining the number of factors in approximate factor models.  Econometrica 70(1), 191–221.
  * Bartlett et al. (2017) Bartlett, P., N. Harvey, C. Liaw, and A. Mehrabian (2017).  Nearly-tight vc-dimension bounds for piecewise linear neural networks.  CoRR abs/1703.02930.
  * Bhadra et al. (2020) Bhadra, A., J. Datta, Y. Li, and N. Polson (2020, August).  Horseshoe Regularisation for Machine Learning in Complex and Deep Models.  International Statistical Review 88(2), 302–320.
  * Brillinger (1977) Brillinger, D. R. (1977).  The identification of a particular nonlinear time series system.  Biometrika 64(3), 509–515.
  * Brillinger (2012) Brillinger, D. R. (2012).  A Generalized Linear Model With “Gaussian” Regressor Variables.  In P. Guttorp and D. Brillinger (Eds.), Selected Works of David Brillinger, Selected Works in Probability and Statistics, pp. 589–606. New York, NY.
  * Carvalho et al. (2012) Carvalho, C. M., H. Lopes, O. Aguilar, and M. Mendoza (2012, 01).  Dynamic stock selection strategies: A structured factor model framework.  Bayesian Stats. 9 9.
  * Chen et al. (2019) Chen, L., M. Pelger, and J. Zhuz (2019, March).  Deep learning in asset pricing.  Technical report, Stanford University.
  * Cheng et al. (2016) Cheng, X., Z. Liao, and F. Schorfheide (2016, 01).  Shrinkage Estimation of High-Dimensional Factor Models with Structural Instabilities.  The Review of Economic Studies 83(4), 1511–1543.
  * Cochrane (2011) Cochrane, J. (2011).  Asset pricing.  Princeton [u.a.]: Princeton Univ. Press.
  * Connor and Korajczyk (1988) Connor, G. and R. Korajczyk (1988).  Risk and return in an equilibrium apt: Application of a new test methodology.  Journal of Financial Economics 21(2), 255–289.
  * Dierckx et al. (2020) Dierckx, T., J. Davis, and W. Schoutens (2020).  Using machine learning and alternative data to predict movements in market risk.
  * Dixon (2022) Dixon, M. (2022).  Industrial forecasting with exponentially smoothed recurrent neural networks.  Technometrics 64(1), 114–124.
  * Dixon et al. (2017) Dixon, M., D. Klabjan, and J. H. Bang (2017).  Classification-based financial markets prediction using deep neural networks.  Algorithmic Finance 6(3-4), 66–77.
  * Dixon and Polson (2020) Dixon, M. and N. Polson (2020).  Short communication: Deep fundamental factor models.  SIAM Journal on Financial Mathematics 11(3), SC–26–SC–37.
  * Dixon and Polson (2022) Dixon, M. and N. Polson (2022).  Bayesian deep fundamental factor models.  In A. Capponi and C.-A. Lehalle (Eds.), Machine Learning in Financial Markets: A guide to contemporary practices. Cambridge: Cambridge University Press.
  * Dobrev and Schaumburg (2013) Dobrev, D. and E. Schaumburg (2013).  Robust forecasting by regularization june 15 , 2013 preliminary and incomplete.
  * Fama and MacBeth (1973) Fama, E. F. and J. D. MacBeth (1973).  Risk, return, and equilibrium: Empirical tests.  Journal of Political Economy 81(3), 607–636.
  * Feng et al. (2018) Feng, G., J. He, and N. G. Polson (2018, Apr).  Deep Learning for Predicting Asset Returns.  arXiv e-prints, arXiv:1804.09314.
  * Frank and Friedman (1993) Frank, I. E. and J. H. Friedman (1993).  A Statistical View of Some Chemometrics Regression Tools.  Technometrics 35(2), 109–135.
  * Golub and Van Loan (2013) Golub, G. H. and C. F. Van Loan (2013).  Matrix computations.  JHU press.
  * Gramacy and Lian (2012) Gramacy, R. B. and H. Lian (2012).  Gaussian process single-index models as emulators for computer experiments.  Technometrics 54(1), 30–41.
  * Gu et al. (2021) Gu, S., B. Kelly, and D. Xiu (2021).  Autoencoder asset pricing models.  Journal of Econometrics 222(1, Part B), 429–450.  Annals Issue:Financial Econometrics in the Age of the Digital Economy.
  * Gu et al. (2018) Gu, S., B. T. Kelly, and D. Xiu (2018).  Empirical asset pricing via machine learning.  Chicago Booth Research Paper 18-04.
  * Hahn et al. (2013) Hahn, P. R., C. M. Carvalho, and S. Mukherjee (2013).  Partial Factor Modeling: Predictor-Dependent Shrinkage for Linear Regression.  Journal of the American Statistical Association 108(503), 999–1008.
  * Harvey et al. (2015) Harvey, C. R., Y. Liu, and H. Zhu (2015, 10).  … and the Cross-Section of Expected Returns.  The Review of Financial Studies 29(1), 5–68.
  * Helland (1990) Helland, I. S. (1990).  Partial Least Squares Regression and Statistical Models.  Scandinavian Journal of Statistics 17(2), 97–114.
  * Higdon et al. (2008) Higdon, D., J. Gattiker, B. Williams, and M. Rightley (2008).  Computer model calibration using high-dimensional output.  Journal of the American Statistical Association 103(482), 570–583.
  * Hoadley (2001) Hoadley, B. (2001).  Statistical Modeling: The Two Cultures: Comment.  Statistical Science 16(3), 220–224.
  * Horel and Giesecke (2019) Horel, E. and K. Giesecke (2019).  Significance tests for neural networks.
  * Hornik et al. (1989) Hornik, K., M. Stinchcombe, and H. White (1989, July).  Multilayer feedforward networks are universal approximators.  Neural Netw. 2(5), 359–366.
  * Kelly and Pruitt (2015) Kelly, B. and S. Pruitt (2015).  The three-pass regression filter: A new approach to forecasting using many predictors.  Journal of Econometrics 186(2), 294–316.  High Dimensional Problems in Econometrics.
  * Kelly et al. (2019) Kelly, B. T., S. Pruitt, and Y. Su (2019).  Characteristics are covariances: A unified model of risk and return.  Journal of Financial Economics 134(3), 501–524.
  * Kingma and Ba (2014) Kingma, D. P. and J. Ba (2014).  Adam: A method for stochastic optimization.
  * Korteweg and Polson (2009) Korteweg, A. and N. G. Polson (2009).  Corporate credit spreads under parameter uncertainty.  Chicago Booth Research Paper Series.
  * Kozak et al. (2018) Kozak, S., S. Nagel, and S. Santosh (2018).  Interpreting factor models.  The Journal of Finance 73(3), 1183–1223.
  * Lettau and Pelger (2020) Lettau, M. and M. Pelger (2020, 03).  Factors That Fit the Time Series and Cross-Section of Stock Returns.  The Review of Financial Studies 33(5), 2274–2325.
  * Liland et al. (2021) Liland, K. H., B.-H. Mevik, and R. Wehrens (2021).  pls: Partial Least Squares and Principal Component Regression.  R package version 2.8-0.
  * Manne (1987) Manne, R. (1987).  Analysis of two partial-least-squares algorithms for multivariate calibration.  Chemometrics and Intelligent Laboratory Systems 2(1), 187–197.
  * Martin and Mahoney (2018) Martin, C. H. and M. W. Mahoney (2018).  Implicit self-regularization in deep neural networks: Evidence from random matrix theory and implications for learning.  CoRR abs/1810.01075.
  * Mhaskar et al. (2016) Mhaskar, H., Q. Liao, and T. A. Poggio (2016).  Learning real and boolean functions: When is deep better than shallow.  CoRR abs/1603.00988.
  * Moritz and Zimmermann (2016) Moritz, B. and T. Zimmermann (2016).  Tree-based conditional portfolio sorts: The relation between past and future stock returns.
  * Naik and Tsai (2000) Naik, P. and C.-L. Tsai (2000).  Partial least squares estimator for single-index models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 62(4), 763–771.
  * Nielsen and Bender (2010) Nielsen, F. and J. Bender (2010).  The fundamentals of fundamental factor models.  Technical Report 24, MSCI Barra Research Paper.
  * Poggio (2016) Poggio, T. (2016).  Deep Learning: Mathematics and Neuroscience.  A Sponsored Supplement to Science Brain-Inspired intelligent robotics: The intersection of robotics and neuroscience, 9–12.
  * Polson and Rockova (2018) Polson, N. and V. Rockova (2018, Mar).  Posterior Concentration for Sparse Deep Learning.  arXiv e-prints, arXiv:1803.09138.
  * Polson et al. (2021) Polson, N., V. Sokolov, and J. Xu (2021).  Deep learning partial least squares.
  * Polson and Scott (2010) Polson, N. G. and J. G. Scott (2010).  Shrink Globally, Act Locally: Sparse Bayesian Regularization and Prediction, Volume 105.  Oxford University Press.  Bayesian Statistics 9.
  * Polson and Scott (2012) Polson, N. G. and J. G. Scott (2012).  Local shrinkage rules, Lévy processes and regularized regression.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 74(2), 287–311.
  * Rosenberg and Marathe (1976) Rosenberg, B. and V. Marathe (1976).  Common factors in security returns: Microeconomic determinants and macroeconomic correlates.  Research Program in Finance Working Papers 44, University of California at Berkeley.
  * Rosipal et al. (2001) Rosipal, R., M. Girolami, L. J. Trejo, and A. Cichocki (2001).  Kernel PCA for Feature Extraction and De-Noising in Nonlinear Regression.  Neural Computing & Applications 10(3), 231–243.
  * Rosipal and Trejo (2001) Rosipal, R. and L. J. Trejo (2001).  Kernel partial least squares regression in reproducing kernel hilbert space.  Journal of machine learning research 2(Dec), 97–123.
  * Simon et al. (2011) Simon, N., J. Friedman, T. Hastie, and R. Tibshirani (2011).  Regularization paths for cox’s proportional hazards model via coordinate descent.  Journal of Statistical Software 39(5), 1–13.
  * Stein (2012) Stein, H. (2012).  Counterparty risk, cva, and basel iii.
  * Stock and Watson (2010) Stock, J. and M. Watson (2010).  Dynamic Factor Models.  Oxford: Oxford University Press.
  * Swanson (2014) Swanson, T. (2014).  Learning from bitcoin’s past to improve its future.
  * Telgarsky (2016) Telgarsky, M. (2016).  Benefits of depth in neural networks.  CoRR abs/1602.04485.
  * Tishby and Zaslavsky (2015) Tishby, N. and N. Zaslavsky (2015).  Deep learning and the information bottleneck principle.  CoRR abs/1503.02406.
  * Tsay (2010) Tsay, R. S. (2010).  Analysis of Financial Time Series (3rd ed.).  Wiley.
  * Wold et al. (1984) Wold, S., A. Ruhe, H. Wold, and W. J. Dunn, III (1984).  The Collinearity Problem in Linear Regression. The Partial Least Squares (PLS) Approach to Generalized Inverses.  SIAM Journal on Scientific and Statistical Computing 5(3), 735–743.
  * Wood (2005) Wood, M. (2005).  Bootstrapped confidence intervals as an approach to statistical inference.  Organizational Research Methods 8(4), 454–470.

##  Appendix A Factor Model Description

Table 3: A short description of the factors used in the Russell $1000$ deep
learning factor model.

ID | Symbol | Value Factors
---|---|---
1 | B/P | Book to Price
2 | CF/P | Cash Flow to Price
3 | E/P | Earning to Price
4 | S/EV | Sales to Enterprise Value (EV). EV is given by
|  | EV=Market Cap + LT Debt + max(ST Debt-Cash,0),
|  | where LT (ST) stands for long (short) term
5 | EB/EV | EBIDTA to EV
6 | FE/P | Forecasted E/P. Forecast Earnings are calculated from Bloomberg earnings consensus estimates data.
|  | For coverage reasons, Bloomberg uses the 1-year and 2-year forward earnings.
17 | DIV | Dividend yield. The exposure to this factor is just the most recently announced annual net dividends
|  | divided by the market price.
|  | Stocks with high dividend yields have high exposures to this factor.
|  | Size Factors
8 | MC | Log (Market Capitalization)
9 | S | Log (Sales)
10 | TA | Log (Total Assets)
|  | Trading Activity Factors
11 | TrA | Trading Activity is a turnover based measure.
|  | Bloomberg focuses on turnover which is trading volume normalized by shares outstanding.
|  | This indirectly controls for the Size effect.
|  | The exponential weighted average (EWMA) of the ratio of shares traded to shares outstanding:
|  | In addition, to mitigate the impacts of those sharp shortlived spikes in trading volume,
|  | Bloomberg winsorizes the data:
|  | first daily trading volume data is compared to the long-term EWMA volume(180 day half-life),
|  | then the data is capped at 3 standard deviations away from the EWMA average.
|  | Earnings Variability Factors
12 | EaV/TA | Earnings Volatility to Total Assets.
|  | Earnings Volatility is measured
|  | over the last 5 years/Median Total Assets over the last 5 years
13 | CFV/TA | Cash Flow Volatility to Total Assets.
|  | Cash Flow Volatility is measured over the last 5 years/Median Total Assets over the last 5 years
14 | SV/TA | Sales Volatility to Total Assets.
|  | Sales Volatility over the last 5 years/Median Total Assets over the last 5 year
|  | Volatility Factors
15 | RV | Rolling Volatility which is the return volatility over the latest 252 trading days
16 | CB | Rolling CAPM Beta which is the regression coefficient
|  | from the rolling window regression of stock returns on local index returns
|  | Growth Factors
7 | TAG | Total Asset Growth is the 5-year average growth in Total Assets
|  | divided by the Average Total Assets over the last 5 years
18 | EG | Earnings Growth is the 5-year average growth in Earnings
|  | divided by the Average Total Assets over the last 5 years
|  | GICS sectorial codes
19-24 | (I)ndustry | $\\{10,20,30,40,50,60,70\\}$
25-35 | (S)ub-(I)ndustry | $\\{10,15,20,25,30,35,40,45,50,60,70,80\\}$
36-45 | (SE)ctor | $\\{10,15,20,25,30,35,40,45,50,55,60\\}$
46-49 | (I)ndustry (G)roup | $\\{10,20,30,40,50\\}$

##  Appendix B Feedforward Network Sensitivities

In a linear regression model

$\hat{Y}=F_{\mathbf{\beta}}(X):=\beta_{0}+\beta_{1}X_{1}+\dots+\beta_{K}X_{K},$
(17)

the model sensitivities are $\partial_{X_{i}}\hat{Y}=\beta_{i}$. In a
feedforward neural network, we can use the chain rule to obtain the model
sensitivities

$J:=\partial_{X}\hat{Y}=W^{(L)}J(I^{(L-1)})=W^{(L)}D(I^{(L-1)})W^{(L-1)}\dots
D(I^{(1)})W^{(1)},$ (18)

where $D_{ii}(I):=\sigma^{\prime}(I_{i}),\leavevmode\nobreak\
D_{ij}=0,\leavevmode\nobreak\ i\neq j$ is a diagonal matrix.

###  B.1 Feedforward network interaction effects

The pairwise interaction effects are readily available by evaluating the
elements of the Hessian matrix. For a L layer network, we define the
$(i,j)^{(th)}$ element of the Hessian of the response,
$\hat{Y}=G_{\hat{\theta}}(X)$, as

$\partial^{2}_{X_{i}X_{j}}\hat{Y}=\sum_{\ell=1}^{L-1}H_{i,j,\ell},\leavevmode\nobreak\
H_{i,j,\ell}:=W^{(L)}D^{(L-1)}W^{(L-1)}\dots\partial_{X_{j}}D^{(\ell)}W^{(\ell)}\dots
w_{i}^{(1)}.$ (19)

where it is assumed that the activation function is at least twice
differentiable everywhere, e.g. $tanh(x)$, softplus etc.

##  Appendix C Proof of Non-linear Stochastic Discount Factors

###### Proof.

Denote the regression coefficients
$a=g(\mathbf{0})\in\mathbb{R},\mathbf{b}=g^{\prime}(\mathbf{0})\in\mathbb{R}^{K},G=g^{\prime\prime}(\mathbf{0})\in\mathbb{R}^{K\times
K}$ are given by the Taylor expansion (to second order):

$g(f)=a+\mathbf{b}^{\top}(\mathbf{f}-\mathbb{E}[\mathbf{f}])+\frac{1}{2}(\mathbf{f}-\mathbb{E}[\mathbf{f}])^{\top}G(\mathbf{f}-\mathbb{E}[\mathbf{f}])+\dots$

$1=\mathbb{E}[mR]=a\mathbb{E}[R]+cov(\mathbf{f}^{\top}R)\Sigma^{-1}\Sigma\mathbf{b}+\frac{1}{2}\mathbb{E}[(\mathbf{f}-\mathbb{E}[\mathbf{f}])^{\top}G(\mathbf{f}-\mathbb{E}[\mathbf{f}])R]+\dots$

Applying the trace trick, then using the invariance of the trace operator
under cyclic permutations, and then the linearity of the trace operator

$\displaystyle\mathbb{E}[(\mathbf{f}-\mathbb{E}[\mathbf{f}])^{\top}G(\mathbf{f}-\mathbb{E}[\mathbf{f}])R]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[tr[(\mathbf{f}-\mathbb{E}[\mathbf{f}])^{\top}G(\mathbf{f}-\mathbb{E}[\mathbf{f}])]R\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[tr[G(\mathbf{f}]-\mathbb{E}[\mathbf{f}])(\mathbf{f}]-\mathbb{E}[\mathbf{f}]^{\top})]R\right]$
$\displaystyle=$ $\displaystyle
tr\left[\mathbb{E}[G(\mathbf{f}-\mathbb{E}[\mathbf{f}])(\mathbf{f}-\mathbb{E}[\mathbf{f}]^{\top})R]\right]$
$\displaystyle=$ $\displaystyle
tr\left[G\mathbb{E}[(\mathbf{f}-\mathbb{E}[\mathbf{f}])(\mathbf{f}-\mathbb{E}[\mathbf{f}]^{\top})R]\right]$

We can then write:

$1=\mathbb{E}[mR]=a\mathbb{E}[R]+cov(\mathbf{f}^{\top}R)\Sigma^{-1}\Sigma\mathbf{b}+\frac{1}{2}\mathbf{b}\Sigma\Sigma^{-1}\mathbf{b}^{-1}tr\left[G\mathbb{E}[(\mathbf{f}-\mathbb{E}[\mathbf{f}])(\mathbf{f}-\mathbb{E}[\mathbf{f}])^{\top}R]\right]\mathbf{b}^{-\top}\Sigma^{-1}\Sigma\mathbf{b}+\dots,$

where $x^{-\top}\equiv[x^{-1}]^{\top}$. This expression can be simplified to

$\mathbb{E}[R]=\alpha+\mathbf{\bm{\beta}}^{\top}\mathbf{\bm{\lambda}}+\frac{1}{2}\mathbf{\bm{\lambda}}^{\top}\Gamma\mathbf{\bm{\lambda}}+\dots,$

where
$\alpha=\frac{1}{a},\mathbf{\bm{\beta}}=cov(\mathbf{f}R)\Sigma^{-1},\mathbf{\bm{\lambda}}=-\Sigma\frac{\mathbf{b}}{a},$
and
$\Gamma=\Sigma^{-1}\mathbf{b}^{-1}tr\left[aG\mathbb{E}[(\mathbf{f}-\mathbb{E}[\mathbf{f}])(\mathbf{f}-\mathbb{E}[\mathbf{f}])^{\top}R]\right]\mathbf{b}^{-\top}\Sigma^{-1}$
and the regression coefficients
$\alpha=f(\mathbf{0})\in\mathbb{R},\mathbf{\bm{\beta}}=h^{\prime}(\mathbf{0})\in\mathbb{R}^{K},\Gamma=h^{\prime\prime}(\mathbf{0})\in\mathbb{R}^{K\times
K}$ are given by the Taylor expansion (to second order):

$g(\mathbf{\bm{\lambda}})=a+\mathbf{\bm{\beta}}^{\top}\mathbb{\bm{\lambda}}+\frac{1}{2}\mathbf{\bm{\lambda}}^{\top}\Gamma\mathbf{\bm{\lambda}}+\dots$

∎

##  Appendix D Additional Results

###  D.1 Factor tilt heatmaps

|
---|---
(a) DPLS(50%) | (b) NN
|
(c) LASSO | (d) Random (WN)
Figure 17: The rescaled factor tilts of equally weighted portfolios
constructed from the predicted top performing 50 stocks against time.

###  D.2 Sector tilt heatmaps

|
---|---
(a) DPLS(50%) | (b) NN
|
(c) LASSO | (d) Random (WN)
Figure 18: The rescaled sector tilts of equally weighted portfolios
constructed from the predicted top performing 50 stocks against time.