# Missing Values and the Dimensionality of Expected Returns

Andrew Y. Chen
Federal Reserve Board Jack McCoy
Columbia University

(May 2023††thanks: First submitted to arXiv.org: July 20, 2022. Emails:
<EMAIL_ADDRESS>(Chen, corresponding author) and
<EMAIL_ADDRESS>(McCoy). This project originated from many
conversations with Fabian Winkler. We thank Heiner Beckmeyer, Charlie Clarke,
Harry Mamaysky, Markus Pelger, Yinan Su, and an anonymous referee for helpful
comments. The views expressed herein are those of the authors and do not
necessarily reflect the position of the Board of Governors of the Federal
Reserve or the Federal Reserve System.)

###### Abstract

Combining many cross-sectional return predictors (for example, in machine
learning) often requires imputing missing values. We compare ad-hoc mean
imputation with several methods including maximum likelihood. Surprisingly,
maximum likelihood and ad-hoc methods lead to similar results. This is because
predictors are largely independent: Correlations cluster near zero and 10
principal components (PCs) span less than 50% of total variance. Independence
implies observed predictors are uninformative about missing predictors, making
ad-hoc methods valid. In PC regression tests, 50 PCs are required to capture
equal-weighted expected returns (30 PCs value-weighted), regardless of the
imputation. We find similar invariance in neural network portfolios.

JEL Classification: G0, G1

Keywords: stock market predictability, stock market anomalies, missing values

##  1 Introduction

A growing literature applies methods from machine learning to asset pricing.
These studies combine the information in dozens, or even hundreds, of cross-
sectional stock return predictors. Freyberger et al. (2020) combine 62; Kelly
et al. (2023) combine 138; and Han et al. (2022) combine up to 193 predictors.
Each of these papers finds that there are economic gains to expanding the set
of predictors beyond the five used in Fama and French (2015).111Lewellen
(2015) finds that linear regressions on 15 predictors, including size, B/M,
momentum, profitability, and investment lead to an annualized Sharpe Ratio of
around 0.8 (excluding micro-caps). Freyberger et al. (2020); Kelly et al.
(2023); and Han et al. (2022) find Sharpe ratios in excess of 1.8, and as high
as 2.4 (excluding micro-caps or value-weighting). Other papers that combine
many predictors include Haugen and Baker (1996); Lewellen (2015); Green et al.
(2017); DeMiguel et al. (2020); Gu et al. (2020); Lopez-Lira and Roussanov
(2020); Kozak et al. (2020); Chen and Velikov (2022); Simon et al. (2022);
Chen et al. (2023).

Buried in this literature is the problem of missing values. When learning from
many predictors, the standard practice of dropping stocks with missing values
is often untenable. For example, applying the standard practice to the 125
most-observed predictors in the Chen and Zimmermann (2022) dataset drops 99.8%
of stocks. So even though imputing missing values may seem dangerous, machine
learning researchers often have no choice but to impute.

We find that simply imputing with cross-sectional averages (as in Kozak, Nagel
and Santosh (2020); Gu, Kelly and Xiu (2020)) does a surprisingly good job of
capturing expected returns. We compare this simple mean imputation with
several alternatives in tests that sort stocks on fitted expected returns (a
la Lewellen (2015)). The expected return models include principal components
regression, gradient-boosted regression trees, and neural networks. For almost
all expected return models, all of the imputations lead to very similar
inferences about mean returns and Sharpe ratios. This invariance holds even
for imputations that employ maximum likelihood via the expectation-
maximization (EM) algorithm.

This invariance comes from the fact that cross-sectional predictors are
largely uncorrelated (Green, Hand and Zhang (2013); McLean and Pontiff (2016);
Chen and Zimmermann (2022)). Almost all cross-sectional correlations lie
between $-0.25$ and $+0.25$, and the first 10 principal components span only
50% of total variance. As a result, observed predictors provide little
information about the missing predictors, and one might as well simply impute
missing values with cross-sectional means. More precisely, our EM imputation
boils down to a set of standardized regression coefficients, and these
coefficients imply that a one standard deviation change in an observed
predictor implies less than a 0.1 standard deviation change in the missing
predictor in the vast majority of cases.

High dimensionality is seen not only in the predictive characteristics, but in
their relationship with future returns. We demonstrate this high
dimensionality in out-of-sample principal components regression tests. Each
month, we form long-short portfolios based on regressions of returns on the
lagged principal components (PCs) of 125 predictors. Using standard principal
component analysis (PCA), dozens of PCs are required to capture the potential
expected returns. Using Huang et al. (2022)’s scaled PCA, which takes on
returns data in the PC estimation, leads to higher mean returns for a given
number of PCs. Nevertheless, scaled PCA still implies a very high
dimensionality: 50 scaled-PCs are required to capture the potential 43%
annualized equal-weighted returns, and 25 scaled-PCs are required to capture
the potential 18% annualized value-weighted returns.

This high dimensionality depends little on whether we use simple mean
imputation or the EM algorithm. Indeed, both methods lead to similar mean
returns for any number of PCs, using either standard or scaled PCA.

A similar invariance is found in a broad array of imputation and forecast
methods. To demonstrate this, we add four imputation methods: (1) industry /
size decile means, (2) last observed, (3) EM on AR1 residuals, and (4) Roweis
(1997)’s practical EM for PCA, and four return forecasts: (1) OLS, (2)
gradient-boosted regression trees, (3) a three-hidden-layer neural network,
and (4) a one-hidden-layer neural network. As in the principal component
regression tests, we form long-short portfolios by sorting stocks on real-time
return forecasts.

With one exception, the forecasts lead to similar mean returns across all
imputation methods. This invariance holds despite the fact that EM on AR1
residuals incorporates time-series information, which Bryzgalova et al. (2023)
argue is important. Indeed, we find a similar invariance when we incorporate
time-series information following Bryzgalova et al. (2023)’s methods. The one
exception to the rule is gradient-boosted regression trees, which appears to
struggle with the degenerate distributions produced by ad-hoc
imputations.222Tree-based algorithms typically offer built-in tree-based
missing value handling, which likely leads to better predictions than imputing
with ad-hoc methods.

Interestingly, ML-style methods frequently _under_ perform OLS. This result
may be surprising given the strong performance of ML in Gu et al. (2020),
Simon et al. (2022), and others. However, it is consistent with the high
dimensionality that we document. High dimensional predictability suggests that
overfitting problems are small, and thus the costs of ML-style regularization
may not exceed the benefits. This cost comes in the form of tuning parameters,
which themselves may be subject to overfitting. Indeed, we find that real-time
tuning of PCR often selects too few PCs, leading to underperformance relative
to OLS. High dimensionality also suggests that the interaction effects found
from tree-based methods may have relatively minor benefits. An important
exception to this underperformance is neural networks, which reliably
outperform OLS by large margins, consistent with the strong performance of
neural networks in Gu et al. (2020) and Simon et al. (2022).

We also provide a detailed description of the origins of missing values in the
Chen and Zimmermann (2022) dataset. The most common driver of missingness is
the requirement of having many years of historical data on the stock. Long-run
reversal is observed for 62% of stocks, and 5-year sales growth is observed
for 40%. These observation rates contrast with size and B/M, which only
require recent data, and are observed for about 98% and 70% of stocks,
respectively. Using more specialized data can drop the observed share below
40% (e.g. R&D data or analyst forecasts), but low observed shares are also
driven by a variety of specialized filters, such as dropping stocks with no
payouts or stocks in particular size quantiles. Thus, missingness is often
driven by the desire to simplify analyses and avoid the modeling of
interactions and non-monotonicity. Modeling these effects would require
digging into the code for the hundreds of individual predictors, an issue we
leave for future research.

In addition to our empirical results, we review the theory around EM
imputation as it relates to cross-sectional asset pricing, including the
confusing “missing at random” (MAR) terminology. Though Rubin (1976) described
MAR as “rather intuitive,” statisticians are still trying to clarify
misinterpretations 50 years later (Little (2021)). We explain how the best
motivation for EM imputation is that it is maximum likelihood assuming
“ignorability” (Rubin (1976)) of the missing data mechanism. Assuming
ignorability is arguably necessary for large scale cross-sectional predictor
data, as the missingness patterns we document are likely too complex to be
modeled convincingly.

We post all of our code at https://github.com/jack-mccoy/missing_data. The EM
imputed data can be found at https://sites.google.com/site/chenandrewy/.

###  1.1 Related Literature

In contemporaneous work, Freyberger, Höppner, Neuhierl and Weber (2022)
(FHNW), Bryzgalova, Lerner, Lettau and Pelger (2023) (BLLP), and Beckmeyer and
Wiedemann (2022) also study missing values in cross-sectional predictor data.
Each paper focuses on a different imputation algorithm, though all avoid
modeling the missingness process. FHNW use moment conditions designed for
cross-sectional predictors, BLLP use Xiong and Pelger (Forthcoming)’s latent
factor method, and Beckmeyer and Wiedemann use a masked language model. Our
paper focuses on EM-style algorithms, which are an important benchmark given
their long history in the study of missing values (Hartley (1958); Orchard and
Woodbury (1972); Rubin (1976); Dempster, Laird and Rubin (1977), etc). Indeed,
many textbooks suggest imputing missing values using EM (Jolliffe (2002);
Little and Rubin (2019); Efron and Hastie (2016)). We also provide a detailed
discussion of the origins of missingness for 125 predictors.

Perhaps most important, we are distinct in providing intuition for why simple
mean imputation seems to perform well in machine learning style applications.
The strong performance of mean imputation is seen in our ML portfolios as well
as FHNW’s regularized non-linear strategies. In contrast, Bryzgalova et al.
(2023) and Beckmeyer and Wiedemann (2022) do not compare their imputations to
simple mean imputation in tests that combine many predictors.

The contemporaneous papers also differ in their predictor and stock return
datasets. In this respect, FHNW is the most similar to ours. We both use 80+
predictors from the CZ open-source dataset and cover all stocks in CRSP.
Beckmeyer and Wiedemann (2022) examine 151 characteristics from the Jensen et
al. (Forthcoming) open source dataset and restrict stock-months to those with
data on at least 20% of the predictors. BLLP examine 45 characteristics from
Kozak et al. (2020) and restrict stocks to those with at least one Compustat
characteristic. Given the size and breadth of these data, and the diversity of
statistical methods, the empirical results in these studies are complementary.

Our results provide a rigorous justification for the ubiquitous use of mean
imputation in machine learning papers. Almost every paper that attempts to
combine large sets of predictors imputes with cross-sectional means or
medians, including Gu, Kelly and Xiu (2020); Freyberger, Neuhierl and Weber
(2020); Kozak, Nagel and Santosh (2020).333Other machine learning style papers
that examine imputing with cross-sectional averages include Green et al.
(2017); DeMiguel et al. (2020); Fulop et al. (2022); Han et al. (2021);
Azevedo et al. (2022); Da et al. (2022); Simon et al. (2022); Zhang and Aaraba
(2022); Jensen et al. (2022). This result is important because the algorithms
in these papers are often complex, so a simple imputation algorithm is helpful
for transparency.

We also add to the literature studying the dimensionality of expected returns.
Weighing stocks equally, Green, Hand and Zhang (2014) and Bessembinder, Burt
and Hrdlicka (2021) find that there are dozens of independent contributors to
expected returns. In contrast, Green, Hand and Zhang (2017) and Lettau and
Pelger (2020) find a much smaller dimensionality in value-weighted analyses.
Our results generally find higher dimensionality than any of these studies,
perhaps because of our use of firm-level predictors and the Chen and
Zimmermann (2022) (CZ) dataset. Firm-level data likely provides more
information about expected returns, and the CZ dataset is significantly larger
than the datasets studied in Green, Hand and Zhang (2017) and Lettau and
Pelger (2020).

Our contribution to these papers is to show that high dimensionality is
little-affected by missing value handling. This invariance is important
because high dimensionality goes against the priors of many economists (for
example, Harvey et al. (2016)). The simplest way to square the empirics with
prior beliefs is to point out a statistical error. Mishandling of missing
values is one potential error, as Green, Hand, and Zhang (2014; 2017) impute
with cross-sectional means, and Bessembinder, Burt and Hrdlicka (2021) build
on the CZ portfolios, which drop missing values. Neither of these methods are
recommended by the statistics or applied literature (Little (1992); Graham
(2009); Little and Rubin (2019)). Our analysis shows that maximum likelihood
has little effect on these conclusions.

These findings complement Lopez-Lira and Roussanov (2020), who document large
Sharpe ratios in strategies that hedge the principal components of individual
stock returns. Taken together, these studies are consistent with the view that
much of cross-sectional predictability arises from exposure to many distinct
sources of information. This “zoo” of information takes time to be impounded
into stock prices (Engelberg, McLean and Pontiff (2018); Bowles, Reed,
Ringgenberg and Thornock (2020)), leading to a zoo of stock return anomalies.

##  2 Missing Values in Cross-Sectional Predictors

We describe the magnitude and origins of missingness among published cross-
sectional predictors. We also explain why imputation is required for large
scale machine learning studies.

###  2.1 Cross-Sectional Predictors Data

Our predictors data begins with the March 2022 release of the CZ dataset. This
dataset contains 207 cross-sectional predictors published in academic
journals.

We drop the 35 discrete predictors because the missing values in this subset
are closely related to the interpretation of these predictors by CZ. Many of
these predictors draw on event studies (e.g. exchange switch), which CZ
interpret as a binary firm characteristic (following McLean and Pontiff
(2016)). Others are based on double-sorts (e.g. analyst recommendation and
short interest). In both cases, CZ use missing values for any stocks that
belong in neither the long or short portfolio, though there are other ways to
interpret these studies.

We drop an additional 15 predictors because they have months with fewer than 2
observations in the 1985-2022 sample. This screen ensures the predictor data
is relatively well-behaved and avoids the use of special techniques to handle
predictor-months with no data.

Last, we keep only the 125 predictors with the most observations in 1985\.
This last screen is required to keep the EM algorithm tractable. This excludes
predictors that focus on very specialized subsets of stocks. For example, it
excludes institutional ownership among high short interest stocks which uses
only stocks in the top 1% of short interest. It also excludes R&D ability,
which uses only stock-months with enough historical R&D data to fit a model of
R&D productivity.

The final list includes not only predictors based on CRSP and Compustat, but
predictors that use IBES analyst forecast data. This inclusion is important,
as analyst forecasts may offer a distinct dimension of predictability.
Unfortunately, it excludes option-based predictors, as option prices are not
well-observed until the 1990s.

###  2.2 Where Does Missingness Come From?

Table LABEL:tab:final-list lists the 125 predictors, ranked by the share of
stocks with observed predictors in 1985, a year that is fairly representative
of the full sample. We explain why these predictors have their corresponding
missingness in what follows.

[Table LABEL:tab:final-list about here]

The most-observed predictors are those that use only the previous month of
CRSP data. These predictors include short-term reversal (streversal),
idiosyncratic volatility (idiorisk), and size, and are observed for at least
98% of stocks.

Requiring a longer history of CRSP data leads to a non-trivial share of
missing values. 12-month momentum is observed for 85% of stocks while medium-
run reversal, which uses the past 18 months of CRSP data, is observed for 80%.
A similar amount of CRSP data is required for coskewness, beta, and two zero
trade predictors (zerotrade and zerotradealt12).

Requiring accounting variables drops the observed share into the 60% to 75%
range. Predictors that fall in this range include the popular B/M,
profitability (roe, gp), and investment (assetgrowth) predictors. Many other
accounting-based predictors fall in this range, including employee growth
(hire), leverage, taxes, external financing (xfin), and accruals.

Rates of observability below 50% are driven by a wide variety of issues. Most
of these predictors require a long history of accounting or CRSP data. The
intangible value predictors (intancfp, intanep, intansp) require 5 years of
past cash flow, earnings, and sales, respectively, and are thus only observed
for about 40% of stocks. Long histories of accounting data are also required
for past 5-years sales growth (meanrankrevgrowth), Titman, Wei and Xie
(2004)’s investment predictor, and composite debt issuance. Return seasonality
in years 11-15 (momseason11yrplus, momoffseasonoff11yrplus) requires 15 years
of past returns, leading to a similar observed share.

But many of the less-observable predictors are missing values simply because
the data items they require are often missing. Predictors that used analyst
forecast data are observed for at most 44% of stocks in 1985, as forecasted
EPS (feps) and analyst revision are only observed for this share. Other
analyst-based predictors have even fewer observations, including analyst
value, and earnings streak. These predictors fall below our 125-most-observed
cutoff, being observed for only 25% and 14% of stocks, respectively. Data
items are also often missing for less-common accounting variables, like
excluded expenses (exclexp), advertising expenditure (adexp), and inventory
growth (invgrowth). Like analyst value and earnings streak, advertising
expenditure and inventory growth are observed for less than 1/3 of firms and
do not make it into Table LABEL:tab:miss-sum.

Still others are missing values due to a diverse range of specialized filters.
Asset tangibility (tang) applies only to manufacturing firms. Payout yield
(payoutyield) drops stocks with non-positive payout. In CZ’s interpretation,
turnover volatility (std_turn) drops stocks in the top two size quintiles, and
operating profitability (operprof) drops stocks in the bottom size tercile.
The original Piotroski F-score (PS) study only examines stocks in the highest
BM quintile and organizational capital keeps only stocks with fiscal-year ends
in December, leading both of these predictors to be dropped from Table
LABEL:tab:miss-sum.

As seen in these details, missingness is sometimes used as a substitute for
modeling non-monotonicity and interaction effects. Modeling these interactions
is highly non-trivial. It requires going into the code behind the individual
predictors, commenting out the missingness assignments, and carefully
introducing new interaction variables if necessary. We leave this work for
future research.

###  2.3 Summarizing Missingness for Individual Predictors

Table LABEL:tab:miss-sum describes missingness at the individual predictor
level. Panel (a) shows the distribution of missing values computed as follows:
For each predictor, we compute the share of stocks that have observed
predictors in selected months between 1970 and 2000\. We then compute order
statistics across predictors.

[Table LABEL:tab:miss-sum about here]

Though many cross-sectional predictability studies begin in 1963, at least 5%
of predictors have no observations until about 1975. A few of these predictors
draw on analyst forecast data, which is not generally observed until 1985
(using IBES).

Missingness reaches a steady state in 1985, as the order statistics are
largely unchanged after this date. In this steady state, the typical predictor
is observed for 70% of stocks, and 75% of predictors are observed for at least
50% of stocks. While this distribution suggests that missing values are a
relatively modest problem, we will see that combining large sets of predictors
compounds the issue.

###  2.4 Summarizing Missingness for Sets of Predictors

When combining predictors, one typically needs to have a numerical value for
an entire set of predictors in each stock-month. This requirement compounds
the magnitudes in Panel (a) of Table LABEL:tab:miss-sum, leading to a
dramatically larger missingness problem.

This problem is illustrated in Figure LABEL:fig:missmap, which shows the
“missingness map” for the most observed 125 predictors in June 1985. Roughly
half of the map is non-shaded, indicating that roughly half of the data is
observed. This does not mean, however, that half of the stocks can be used in
a machine learning algorithm without imputation. In general, only stocks with
no missing values can be used, corresponding to columns in the missingness map
that are entirely non-shaded. Almost no columns satisfy this requirement, and
thus one often _must_ impute data in order to apply machine learning in a
meaningful way.

[Figure LABEL:fig:missmap about here]

Panel (b) of Table LABEL:tab:miss-sum quantifies this problem. This panel
examines the prevalence of missing values when combining $J$ predictors, where
$J$ ranges from 25 to 125. The $J=25$ row examines the 25 predictors with the
most observations a given month. It shows that one can combine 25 predictors
with relatively little problem because they are so well-observed. In June
1985, 93.1% of predictor-stocks have observations and 77.5% of stocks have
predictor data for all 25 predictors. Thus, by simply dropping missing values,
one can likely provide a representative description of the stock market for
these 25 predictors.

But when combining for 75 or more predictors, one needs to take a stand on
imputation. As seen in the $J=75$ row, simply dropping stocks with missing
values would result in dropping all but 15.7% of the data. Indeed, for
$J=125$, essentially the entire dataset is gone if one follows the traditional
missing value handling used in individual predictor studies.

Returning to Figure LABEL:fig:missmap, the complex missingness patterns
illustrate why an algorithm like EM is required. There are dozens of distinct
shapes in Figure LABEL:fig:missmap, suggesting that dozens of equations are
required to model the missingness mechanism. It is unlikely that an analyst
could convincingly present such a complicated model. The only alternative is
to assume Rubin ignorability, as in EM imputation.

##  3 EM Imputation

The missingness patterns imply that one _must_ impute data when combining more
than 100 predictors. These patterns are complex and require a general
algorithm that handles arbitrary patterns. EM is one such general algorithm.

In short, our EM imputation assumes that each cross-section of predictors is
drawn from a multivariate normal distribution, and then uses maximum
likelihood to estimate the expected value of the missing predictors. This
section details our implementation of EM (Section 3.1-3.3) and clarifies
common misconceptions about the “missing at random” assumption (Section 3.4).

###  3.1 Missing Data Imputation Under “Ignorability”

At its core, EM imputation assumes that that the missing value mechanism is
“ignorable” in the sense of Rubin (1976). To formalize this idea, suppose that
in month $t$, stock $i$ has a vector of predictors $X_{i,t}$ drawn from a
distribution

$\displaystyle X_{i,t}$ $\displaystyle\sim f\left(x_{i,t}|\theta_{t}\right)$
(1)

where $x_{i,t}$ is a vector of predictor realizations and $\theta_{t}$ is a
vector of parameters. This form means that each month of data is modeled
separately. It also means that a stock’s observed predictors are potentially
informative about its unobserved predictors (high B/M could mean high tail
risk) and is similar to the third moment condition in Freyberger et al.
(2022). For ease of exposition, we assume independence across $i$, but this
assumption can likely be weakened to allow for weak dependence as in
Wooldridge (1994).

EM allows for general missingness patterns, which results in somewhat unusual
notation. Let $x_{i,t}^{obs}$ be a vector that replaces missing values in
$x_{i,t}$ with a non-numeric value $\ast$. Similarly, $x_{i,t}^{miss}$
replaces observed values in $x_{i,t}$ with $\ast$. Define the random vector
counterparts $X_{i,t}^{miss}$ and $X_{i,t}^{obs}$ similarly.

This notation allows us to write down a quasi-likelihood function given only
the observed data

$\displaystyle
L_{ign}\left(\theta_{t}|\left[x_{1,t}^{obs},...,x_{N_{t},t}^{obs}\right]\right)$
$\displaystyle\equiv\Pi_{i=1}^{N_{t}}\int
dx_{i,t}^{miss}f\left(\left[x_{i,t}^{obs},x_{i,t}^{miss}\right]|\theta_{t}\right).$
(2)

where $\int dx_{i,t}^{miss}$ indicates integrating over the support for each
of stock $i$’s missing predictors. The fundamental assumption of EM imputation
is that maximizing Equation (2) provides valid inference about $\theta_{t}$.

The subscript “ign” refers to Rubin (1976) ignorability. Equation (2) is not
the full likelihood function, as it ignores the process that determines
whether, for a particular stock $i$, a particular element of $X_{i,t}$ is
missing. As shown in Section 2.2, this process is driven by several distinct
mechanisms and the resulting missingness patterns are complex. As such,
ignorability is likely a necessary assumption.

Even assuming ignorability, maximizing Equation (2) is not straightforward.
This optimization can be achieved with the EM algorithm of Dempster, Laird and
Rubin (1977) (see also Hartley (1958); Orchard and Woodbury (1972)). EM
iterates between the “expectation” step, which computes a kind of expected
quasi-likelihood assuming some parameter values, and the “maximization” step,
which chooses new parameter values by maximizing the expected quasi-likelihood
from the expectation step. We provide further details on the EM algorithm in
Appendix A.2.

Given $\hat{\theta}_{t}$ that maximizes Equation (2), we impute missing values
with their expected values conditional on observed predictors

$\displaystyle\hat{x}_{i,t}^{miss}$
$\displaystyle=E\left(X_{i,t}^{miss}|X_{i,t}^{obs}=x_{i,t}^{obs},\hat{\theta}_{t}\right).$
(3)

Since the EM algorithm outputs Equation (3) as a byproduct of the expectation
step, it serves as a natural method for imputation.

###  3.2 Multivariate Normal Implementation

The EM algorithm can be quite difficult to implement. In principle, each
expectation step requires an extremely large number of high-dimensional
integrals, and each maximization step requires solving a global optimization
problem.

To keep the problem tractable, we assume that predictors are multivariate
normal after a transformation. That is, we transform each predictor by
winsorizing symmetrically at the 1% level and then apply Hawkins and Weisberg
(2017)’s extension of the Box and Cox (1964) method for making each predictor
approximately normal. Then we standardize predictors to have zero mean and
unit variance. Finally, we assume

$\displaystyle f\left(x_{i,t}|\theta_{t}\right)$
$\displaystyle\sim\text{MVN$\left(0,\Sigma_{t}\right)$}$ (4)

where $\text{MVN}\left(0,\Sigma_{t}\right)$ is a multivariate normal
distribution with mean zero and covariance matrix $\Sigma_{t}$ and $x_{i,t}$
represents the transformed predictors.

While the MVN assumption may appear restrictive, it can be thought of as an
approximation of a more general distribution, and simulation studies often
find that the MVN assumption works well even if the underlying data is not
normal (Little (1988); Azen et al. (1989); Graham and Schafer (1999)). Little
(1988) shows that allowing for fat tails leads to the same algorithm but with
smaller weights on the extreme observations. As we winsorize our data before
imputation, such a fat tail extension should have little effect on our
results.

The Box-Cox transformation is not commonly found in asset pricing, but it is
just an extension of the log transformation that is commonly applied to size
and B/M (Fama and French (1992)). More broadly, cross-sectional predictability
studies rarely impose a functional form on the relationship between the
predictor predictors and expected returns, so one might as well transform the
predictors to be approximately normal.

Importantly, our implementation of EM avoids look-ahead bias as each month of
data is handled independently. Further details on the implementation are found
in Appendix A.1 and A.2.

###  3.3 Predictor Selection for EM

To make the algorithm tractable, we limit our EM analysis to the 125 most-
observed predictors in 1985. (Table LABEL:tab:final-list).

This restriction is important because, even with the MVN assumption, the EM
algorithm is still computationally intensive due to complex missingness
patterns. While the algorithm boils down to a sequence of matrix operations,
every missingness pattern seen in Figure LABEL:fig:missmap implies a different
matrix inversion, and every matrix inversion must be done in every iteration
of the algorithm. As a result, it can take up to four hours to impute a single
month of data for the 125 most-observed predictors. This algorithm must then
be repeated for every month of the data. We make this computation feasible by
imputing all 432 months of data in parallel on a cluster.

Table LABEL:tab:final-list shows that the 125 most-observed predictors span
most economic and data categories. Importantly, our list uses data from not
only CRSP and Compustat, but also IBES. Including IBES data is important, as
it may offer a distinct dimension of predictability not seen in CRSP or
Compustat. The predictors omitted from our list are largely those that only
apply in specialized subsets of firms, like the Piotroski F-score (see Section
2.2).

###  3.4 Discussion of “Missing at Random” and Related Terms

It is well-known that “missing at random” (MAR) is a sufficient condition for
EM imputation (Little and Rubin (2019)). MAR and related terminology
frequently arises in missing data analysis (Little and Rubin (2019); Xiong and
Pelger (Forthcoming)). We briefly discuss MAR and how it relates to our
analysis here.

An intuitive definition of MAR is that whether or not a stock-month is missing
is random:

$\displaystyle
p_{MISS_{i,t}|X_{i,t}}\left(miss_{i,t}|\left\\{x_{i,t}^{obs},x_{i,t}^{miss}\right\\},\phi_{t}\right)$
$\displaystyle=\pi\quad\text{for all }x_{i,t}^{obs},x_{i,t}^{miss},\phi_{t}$
(5)

where $MISS_{i,t}$ is a random vector indicating the missing values for stock
$i$ in month $t$, $miss_{i,t}$ is its realization, $\phi_{t}$ is a parameter
vector that governs the missing value process, and
$p_{MISS_{i,t}|X_{t}}\left(\cdot\right)$ is the probability of missingness
given information on all realized predictors for stock $i$. Indeed, Equation
(5) is how Xiong and Pelger (Forthcoming) defines MAR.

Equation (5), however, is not the definition used in much of the missing data
literature. Graham (2009), Little (2021), and others describe Equation (5) as
“missing _completely_ at random.” These papers follow Rubin (1976) and define
MAR as

$\displaystyle
p_{MISS_{i,t}|X_{i,t,}}\left(miss_{i,t}|\left\\{x_{i,t}^{obs},x_{i,t}^{miss}\right\\},\phi_{t}\right)$
$\displaystyle=p_{MISS_{i,t}|X_{i,t}}\left(miss_{i,t}|\left\\{x_{i,t}^{obs},\tilde{x}_{i,t}^{miss}\right\\},\phi_{t}\right)$
(6) for all $x_{i,t}^{miss}$, $\tilde{x}_{i,t}^{miss}$, and $\phi_{t}$,

where $\tilde{x}_{i,t}^{miss}$ is a realization of the missing predictors. As
seen in Equation (6), Rubin’s MAR does not __ mean that missingness is random
in the sense that it is determined by coin flips. “Missing not at random” is,
then, when Equation (6) is not satisfied.

What does Rubin’s MAR mean? Perhaps the best description is that it is “the
weakest general condition under which pure likelihood inference can be based
on [Equation (2)] without needing a model for the missingness mechanism”
(Little (2021)). This result is important because Equation (6) is a much
weaker assumption than Equation (5). Thus, the EM algorithm does not require
that missingness is completely random. Moreover, EM does allow for missingness
to depend on the data.

##  4 Main Predictability Results

Having described the need for EM imputation, we are finally in a place to
present our main results. This section compares the performance of EM vs mean
imputation for predicting long-short returns.

We begin with a warm-up exercise that studies portfolios sorted on individual
predictors (Section 4.1). Section 4.2 presents the main results, which
examines long-short portfolios that the principal components of 125
predictors.

Though we focus on two imputations and two forecasting methods, Section 5
shows that our main results can be extended to many more methods, including
imputations that use time-series information and forecasts based on neural
networks.

###  4.1 Single-Predictor Strategies

We start with single-predictor strategies because they are so familiar to
asset pricing researchers. It is unlikely that missing data handling has a
large effect here, but these results are helpful “descriptive” statistics
before the main results.

Our single-predictor strategies go long the 500 stocks with the strongest
predictor and short the 500 stocks with the weakest predictor each month. We
do this for two missing data handling methods: (1) EM imputation and (2)
dropping stocks with missing values.

This design allows us to examine the most common method for handling missing
values (dropping missing) while simultaneously removing “predictor dilution”
effects that can occur when using imputed data. Intuitively, the imputed data
covers more stocks, and if these additional stocks have more neutral
predictors, the extreme deciles using imputed data will display less
predictability simply because of this additional coverage. Fixing the number
of stocks in each leg controls for this issue. For simplicity, we drop
predictor-months that have less than 1000 stock-predictor observations. This
filter removes 1% of predictor-months.

Table LABEL:tab:one-signal shows the resulting distribution of mean returns
and Sharpe ratios. EM generally improves performance, on average. However, the
effects of EM imputation are generally small. In equal-weighted portfolios, EM
increases the median annualized return from 5.9% to 7.0% and the median
annualized Sharpe Ratio from 0.53 to 0.56. In value-weighted portfolios, the
median performance is almost unchanged. EM increases somewhat the dispersion
of performance measures, as seen in the slightly larger right and left tails.

[Table LABEL:tab:one-signal about here]

The “EM Improvement” row computes the difference between EM and missing
dropped performance at the predictor level. The typical difference is tiny, at
around 0.1 percentage point per year for expected returns. The largest
difference is the 90th percentile for EM improvement for equal-weighted long
short returns, and this difference is a rather moderate 3.6 percentage points
per year.

###  4.2 Strategies that Combine the Principal Components of 125 Predictors

Here we examine the most compelling setting for studying missing values:
strategies that combine information from many predictors. We combine
information using two kinds of principal components regression (PCR). PCR is a
simple method for handling correlated predictors and potential overfitting
problems. It also provides a natural connection with the dimensionality of
expected returns.

Our PCR analysis forms “real-time” long-short strategies as follows: For each
month $t$ between January 1995 and December 2020, we

  1. 1.

Find the principal components (PCs) of the predictors using data from the past
120 months. We find PCs two ways:

    1. (a)

Standard PCA: find the pooled covariance matrix of $X_{i,t}$ and compute the
$k$th PC as the eigenvector of the covariance matrix with the $k$th largest
eigenvalue.444Pooling can be justified if
$\left\\{N^{-1}\sum_{i=1}^{N}X_{i,t}X_{i,t}\right\\}_{t=1,2,3,...}$ is weakly
dependent, so that the pooled estimator converges to
$E\left(X_{i,1}X_{i,1}^{\prime}\right)$.

    2. (b)

Scaled PCA (Huang et al. (2022)): For each predictor $j$, estimate
$r_{i,t}=\alpha_{j}+\gamma_{j}X_{i,j,t-1}+\varepsilon_{i,j,t}$ using OLS, and
define the scaled predictor $\dot{X}_{i,j,t}\equiv\hat{\gamma}_{j}X_{i,j,t}$.
Then stack into vectors
$\dot{X}_{i,t}\equiv\left[\dot{X}_{i,1,t},...,\dot{X}_{i,J,t}\right]^{\prime}$
and apply standard PCA to get scaled-PCs.

  2. 2.

Regress stock-level returns on lagged first $K$ PCs using OLS. We examine $K$
ranging from 1 to 75.

  3. 3.

Go long the top decile of predicted returns and short the bottom decile.

  4. 4.

Hold the position for one month.

We apply this procedure using either EM imputation or simple imputation with
cross-sectional means. As seen in Section 2.4, _some_ form of imputation is
required when combining 125 predictors.

####  4.2.1 Regressions on Standard PCs

Figure LABEL:fig:pcr shows PCR results using standard PCs. It plots the mean
returns (Panel (a)) and Sharpe ratios (Panel (b)) of strategies as a function
of the number of PCs used to forecast returns. The various lines show equal-
or value- weighting, and EM or simple mean imputation.

[Figure LABEL:fig:pcr about here]

Notably, the results depend very little on whether EM or mean imputation is
used. Echoing Table LABEL:tab:one-signal, EM imputation (solid line) leads to
slightly better performance than simple mean imputation (dashed line), but the
improvement is rather modest. The similarity in performance is seen throughout
the plot, for all choices for the number of PCs. This result suggests that the
simple mean imputation used in the machine learning strategies of Freyberger
et al. (2020); Gu et al. (2020) and others does not miss out on the potential
returns. We verify this robustness for some machine learning strategies in
Section 5.

Also notably, mean returns and Sharpe ratios increase in the number of PCs
throughout the plot, all the way up to 60 PCs. The improvements are large. The
mean return using 60 PCs is nearly twice as large as the mean returns from
using 20 PCs using equal- or value-weighting. Similar magnitudes are seen in
the Sharpe ratios.555In a previous version of this paper, we found that using
“value-weighted” OLS in PCR results in about 30 PCs being required to capture
the potential value-weighted returns. To our knowledge, value-weighted OLS
lacks a solid econometric foundation, so we opted to simplify and focus on
standard OLS in this draft

This result points to a much higher dimensionality of expected returns than is
commonly believed among asset pricing researchers. Going back at least to Fama
and French (1993), the asset pricing literature has characterized the cross-
section of expected returns as being summarized by just a few factors. More
recently these three factors have been expanded to include five or six total
factors (Fama and French (2015)). But this number is still far fewer than the
60 dimensions of predictability seen in Figure LABEL:fig:pcr.

####  4.2.2 Regressions on Scaled PCs

A potential explanation for this apparent high dimensionality is that standard
PCA does not incorporate target information when extracting PCs. Huang et al.
(2022)’s scaled-PCA provides an intuitive method for incorporating this
information. They show that simply scaling predictor variables by their
covariance with the target variable leads to PCR forecasts that outperform
standard PCR. They demonstrate this outperformance using theory, simulation,
and forecasting using macroeconomic data.

Figure LABEL:fig:spca shows the result of using scaled-PCs in our PCR tests.
Like Huang et al. (2022), we find that scaled-PCA is more efficient than PCA.
The improvement is large: 10 scaled PCs leads to equal-weighted returns of
about 30% per year, nearly double the roughly 18% found using standard PCs.
Outperformance is also seen in value-weighted portfolios, where 10 scaled PCs
leads to returns of around 11%, while 10 standard PCs leads to returns of
around 8%.

[Figure LABEL:fig:spca about here]

While scaled-PCA implies that fewer PCs are required, it still implies a very
high dimensionality. At least 50 scaled PCs are required to capture the
potential equal-weighted returns, while at least 25 scaled PCs are required to
capture the potential value-weighted returns. As with standard PCA, this
inference depends little on the missing data handling.

###  4.3 Imputation and Dimensionality

The strong performance of ad-hoc imputation is closely related to
dimensionality. To see the intuition, note that under the multivariate normal
assumption (Equation (4)), imputation (Equation (3)) boils down to a
regression equation

$\displaystyle\hat{x}_{i,j,t}^{miss}$
$\displaystyle=\hat{\beta}_{i,j,t}^{\prime}x_{i,t}^{obs}$ (7)

where $x_{i,j,t}^{miss}$ is a component of the vector $x_{i,t}^{miss}$ and the
imputation slope vector $\hat{\beta}_{i,j,t}^{\prime}$ is given by

$\displaystyle\hat{\beta}_{i,j,t}^{\prime}$
$\displaystyle=\hat{\Sigma}_{j,obs|i,t}\hat{\Sigma}_{obs|i,obs|i,t}^{-1}$ (8)

where $\hat{\Sigma}_{obs|i,obs|i,t}$ is the submatrix of $\hat{\Sigma}_{t}$
corresponding to the observed predictors for stock $i$, and
$\hat{\Sigma}_{j,obs|i,t}$ is the subvector of $\hat{\Sigma}_{t}$
corresponding to the covariance between predictor $j$ and the observed
predictors.

In essence, our imputation forecasts missing values (LHS of Equation (7))
based on their covariance with observed values (RHS of Equation (8)), just as
in OLS. If there was no missing data, Equation (8) would just be the slope
from regressing predictor $j$ on all of the other predictors. However, missing
data present a non-trivial problem for OLS, as it is not at all clear how one
should compute covariances when arbitrary subsets of the data are missing. The
purpose of EM and the associated formalism is to handle this problem. But in
the end, the intuition from OLS still holds. If the LHS variable is
uncorrelated with the RHS variable, then the slope of their relationship is
zero. In this extreme case, imputing with cross-sectional means is no
different than OLS, and indeed may be more efficient.

More generally, imputation will be similar to simple mean imputation if the
off-diagonal terms of $\Sigma_{t}$ cluster close to zero. Similarly, if
$\Sigma_{t}$ is high-dimensional (the eigenvalues of $\Sigma_{t}$ are similar
in magnitude), then the predictors are distinct and provide limited
information about each other. In either case, intuition suggests that EM
imputation should have little effect. Deriving a formal relationship between
the distribution of off-diagonal terms and the eigenvalues of $\Sigma_{t}$ is
difficult, but we provide a simple simulation that illustrates this
relationship in Appendix A.3.

Figure LABEL:fig:cor examines the off-diagonal terms of $\Sigma_{t}$ in the
empirical predictor data. The “Observed” lines simply calculate correlations
between two predictors using all pairs of stocks that have predictor data in
the given month (sometimes called “available case” or “pairwise complete”
correlations). The left panels show that these observed correlations cluster
close to zero, consistent with Chen and Zimmermann (2022). This distribution
suggests a very large number of PCs are required to describe the data, and
this result is verified in the right panels. Based on observed correlations,
more than 9 PCs are required to capture just half of total variance.

[Figure LABEL:fig:cor about here]

Available case analysis is not generally a recommended solution from the
missing data literature (e.g. Little (1992)). The “EM Algo” lines in Figure
LABEL:fig:cor examine a recommended solution: correlations estimated from the
EM algorithm. Despite the additional rigor, Figure LABEL:fig:cor shows that EM
leads to a nearly identical distribution (Appendix A.4 shows that the
distributions are not exactly identical). The similarity is reminiscent of
Xiong and Pelger (Forthcoming), who find that available case covariance
matricies are consistent estimators under some assumptions on the missing data
process. This result also makes sense given that the most significant problems
with available case analysis seem to occur when the data are highly correlated
(Haitovsky (1968)).

The panels in Figure LABEL:fig:cor show that high dimensionality is stable
over time. Regardless of whether one computes the correlations in June 1990,
June 2000, or June 2020, almost all correlations lie between -0.25 and +0.25,
10 PCs span only 50% of total variance, and both the EM algorithm and observed
correlations lead to nearly identical distributions and PCA decompositions.

Given this high dimensionality, Equation (8) implies that imputation slopes
will be close to zero. This idea is confirmed in Figure LABEL:fig:slopes,
which shows the distribution of imputation slopes in selected months between
1990 and 2010. The slopes are computed by plugging in the EM-estimated
covariance matrix into Equation (8), and then stacking all of the slope
vectors into one long vector before plotting the distribution.

[Figure LABEL:fig:slopes about here]

Figure LABEL:fig:slopes shows that almost all slopes lie between -0.2 and
+0.2, regardless of the month one examines. Indeed the mean absolute slope is
just 0.065. As our predictor data is standardized to have unit variance, these
coefficients imply that a one standard deviation increase in an observed
predictor implies a 6.5% change to the missing predictor—almost no change at
all. As a result, simple mean imputation and EM imputation will lead to
similar results for examining expected returns.

##  5 Comparing 6 Imputations using 6 Forecasting Methods

The thrust of this paper examines just two imputation methods and two return
forecasting methods. This section expands the analysis in both directions.

We examine four additional imputations, including ones that incorporate time-
series information. We examine four additional return forecasting methods,
including two based on neural networks.

The main results continue to hold: ad-hoc imputations and more sophisticated
imputations lead to largely similar inferences about expected returns.

###  5.1 Four Additional Imputation Methods

We examine the following additional imputations:

  1. 1.

MVN EM on AR1 Residuals: We estimate an AR1 model on Box-Cox transformed
predictors using the past 60 months of data and apply MVN EM to the residuals.
To impute missing predictors, we add the AR1 prediction to the MVN-EM-imputed
residuals. This algorithm is consistent (Appendix A.5.1) and takes on time-
series information, which BLLP find is important in their imputation error
tests.

  2. 2.

Practical EM for Probabilistic PCA: We assume Box-Cox transformed predictors
follow a 60-dimensional factor structure and impute using Roweis (1997)’s
“practical” EM algorithm for missing data. This practical EM algorithm simply
imputes missing predictors using the factor model prediction in the E-step,
and as such is not a maximum likelihood estimate. However, it is easily
implemented using the pcaMethods R library and provides a regularized
alternative to our MVN EM imputation.

  3. 3.

Industry - Size Decile Means: Within each 3-digit SIC code, we sort stocks
into size deciles. We then replace missing Box-Cox transformed predictors with
their means within each industry-size group. This imputation takes on the
intuition that industry and size are important determinants of a firm’s
economics.

  4. 4.

Last Observed: Replace missing values with the last observed value within the
last 12 months. If there is no observed value within the last 12 months, we
replace missing values with cross-sectional means. This imputation takes on
time-series information while avoiding using stale information. It is also
easy to implement

For additional details see Appendix A.5.

###  5.2 Four Additional Return Forecasting Methods

We examine the following additional return forecasting methods:

  1. 1.

Simple OLS: We simply regress future stock returns on predictors using OLS
without any regularization. Simple OLS provides a helpful benchmark.

  2. 2.

Gradient-Boosted Regression Trees (GBRT): We model future stock returns using
as the average of $B$ shallow regression trees. The model is fitted by
gradient boosting: starting from an existing tree, the tree is “boosted” by
adding a new tree, where the new tree is chosen to step along the gradient of
a squared-error loss function, and this process is continued $B$ times. $B$
and other hyperparameters are tuned following Gu et al. (2020) and Li et al.
(2022), and the algorithm is implemented with Microsoft’s open source lightGBM
library.

  3. 3.

3-Layer Neural Network (NN3): We model the relationship between future stock
returns and predictors as the average of an ensemble of feed-forward neural
networks with three hidden layers. We fit the model using an Adam optimizer
with batch normalization and early stopping. Tuning and other specification
details follow Gu et al. (2020), who find that neural network performance
peaks at three hidden layers. We implement the algorithm with Google’s open
source Tensorflow and Keras libraries.

  4. 4.

1-Layer Neural Network (NN1): This forecast uses the same model as NN3, but
with only one hidden layer. This model highlights the robustness of NN3, which
we find is by far the best-performing forecasting method.

For additional details on the neural network specifications see Appendix A.6.

###  5.3 Out-of-Sample Forecasting and Hyperparameter Tuning

We impute missing values using one of the six methods and then construct “out-
of-sample” portfolios using one of the six forecasting methods.

Portfolio construction uses the following procedure. For each June $t$ between
1995-2019, we

  1. 1.

Separate the data between January 1985 and $t$ into a training sample and a
validation sample. The validation sample is either the latter half of the data
or the past 12 years, whichever sample is shorter. The remainder of the data
is the training sample.

  2. 2.

For all hyperparameter sets in a pre-selected set of sets, fit returns to
lagged imputed data using the training sample and measure the RMSE in the
validation sample.

  3. 3.

Using the hyperparameter set with the smallest RMSE, fit returns to lagged
imputed data from January 1985 to $t$.

  4. 4.

For each month $t,t+1,...,t+11$, sort stocks on the predicted return and form
long-short decile portfolios.

This procedure follows Gu et al. (2020) as closely as possible. We differ only
in requiring that at least half of the January 1985 to $t$ subsamples are used
for training. Gu et al. do not need this requirement because their data goes
back to 1957.

With the exception of OLS, the forecasting methods require hyperparameter
tuning. We tune hyperparameters following Gu et al. (2020) as closely as
possible (see their Appendix Table A.5). For PCR and sPCR, we select the
number of PCs from $\left\\{10,30,50,70,90\right\\}$. For GBRT, the depth is 1
or 2, the learning rate is 0.01 or 0.10, and the number of trees is in
$\\{1,250,500,750,1000\\}$. For NN1 and NN3, the L1 penalty is either
$10^{-5}$ or $10^{-3}$, the learning rate is 0.001 or 0.010, and other
hyperparameters are fixed: batch size is 10,000, epochs is 100, early stopping
patience is 5, the ensemble number is 10, and we use the Adam optimizer with
default parameters in the Keras interface to TensorFlow.666We differ from Gu
et al. (2020) by using the standard L2 objective in GBRT rather than a Huber
loss. Our implementation of GBRT uses lightGBM, which currently has poor
documentation of the objective function, and we find that objective = ’huber’
works poorly on simulated and empirical data. Li et al. (2022) does not
mention using Huber loss.

###  5.4 Results Across a Broad Array of Methods

Table LABEL:tab:imp-fore-grid shows the result across all combinations of the
6 imputations (rows) and 6 forecasting methods (columns). We examine both
equal-weighted (Panel (a)) and value-weighted (Panel (b)) returns.

[Table LABEL:tab:imp-fore-grid about here]

####  5.4.1 Invariance Across Imputations

The central feature of the table is that mean returns vary little across
imputations. This invariance is seen in the fact that returns change very
little when one moves down the rows of the table. For example, in Panel (a)
(equal-weighted returns) the OLS column ranges from 47% (using ad-hoc
imputations) to 49% (using EM-imputations). Panel (b) shows a similar result
for OLS under value-weighting, where returns vary from 20% to 22% (using ad-
hoc methods) to 19% to 25% using EM.

EM-style methods generally lead to higher returns but the improvements are
modest and somewhat uncertain. Using OLS, PCR, sPCR, NN3, or NN1, EM-style
methods consistently outperform by a couple percentage points under equal
weighting, but the outperformance is variable under value-weighting. For
roughly half of the EM-based imputations, value-weighted returns are slightly
smaller than those obtained from simple mean imputation. Regardless, the main
message of the table echoes Figures LABEL:fig:pcr and LABEL:fig:spca: formal
imputation has little impact on inferences about mean returns.

Notably, imputing with MVN on AR1 residuals leads to similar results as simple
mean imputation. This finding illustrates the robustness of our main results
to the incorporation of time-series information, which BLLP find is important
for imputation accuracy under simulated missingness patterns. We also find
that the BLLP local backward cross-sectional imputation has little effect on
estimated expected returns (Appendix Figure A.3). We view this result as
intuitive since missingness tends to occur in large blocks in the time
dimension. If a stock is missing B/M in a particular month, we find that is
also missing B/M 12 months ago with 93% probability, and that it is missing
B/M 24 months ago with 96% probability. We find similar time-series blocks of
missingness for momentum and asset growth. Indeed, MVN on AR1 residuals leads
to slightly lower mean returns than pure cross-sectional EM, perhaps because
the amount of time-series information is too small to overcome the fitting
error that comes from imposing a specific form of AR1 model.

An exception to this invariance is GBRT. Using GBRT, EM-style imputation leads
to noticeably higher returns under equal- or value-weighting, perhaps
regression trees deal poorly with the degenerate distributions that result
from ad-hoc imputation. One should not over-interpret this result, however,
since tree-based algorithms often have their own specialized methods for
handling missing values, which we do not employ because of our goal of
evaluating “general-purpose” imputations.777lightGBM can handle missing values
by adding a second greedy stage to each split assessment. This mechanism is
not well-documented but is described on Github
(https://github.com/microsoft/LightGBM/issues/2921). For random forests,
Brieman suggests using a previously fitted tree to impute missing values based
on “proximities” with the non-missing value cases
(https://www.stat.berkeley.edu/$\sim$breiman/RandomForests/cc_home.htm#missing1).

####  5.4.2 Performance Variation Across Forecasting Methods

A second theme from Table LABEL:tab:imp-fore-grid is that ML-style methods
often _under_ perform relative to simple OLS. Moving across the columns of
Table LABEL:tab:imp-fore-grid, outperformance relative to OLS is only seen
using neural network forecasts.

These results may be surprising, given that many papers tout the
outperformance of machine learning methods (Gu et al. (2020); Freyberger et
al. (2020); Simon et al. (2022)). However, they are consistent with these
papers once one examines the details. Gu et al. (2020) find that applying PCR
and GBRT to 94 predictors leads to returns that are not much larger than those
obtained from applying linear regression to just three predictors. Freyberger
et al. (2020) and Simon et al. (2022) both examine continuous non-linear
models of returns, which are differ substantially from linear PCR models and
the discontinuous GBRT models.

Moreover, the underperformance of PC-based methods is a natural consequence of
the high dimensionality of expected returns (Section 4.2). Unlike OLS, which
requires no tuning, PC-based methods requires the choice of the number of PCs
to use. Since predictability is high dimensional (Figures LABEL:fig:pcr and
LABEL:fig:spca), the optimal choice of PCs is very large, at 25 or more for
value-weighted portfolios, and 50 or more for equal-weighted portfolios. These
optimal choices are not known in real-time, however. To analyze real-time
returns, the standard practice is tune the number of PCs to minimize errors in
a rolling validation sample (Gu et al. (2020)). We find that this kind of
tuning often leads to choosing only 10 PCs, and thus suboptimal returns.
Intuitively, high dimensionality of returns implies that the cost of
overfitting is small, and thus regularization may be suboptimal.

A final theme from Table LABEL:tab:imp-fore-grid is that neural network models
reliably outperform other methods. Mean returns are at least 20 percent larger
than other methods using equal-weighting, and about 50 percent larger using
value-weighting. In fact, we find that other neural network architectures work
quite well also, and that performance is fairly robust to ensembling and other
tuning choices if early stopping is set following Gu et al. (2020). These
results support the strong performance of neural networks found in other
papers (Gu et al. (2020); Simon et al. (2022)) and shows that this strong
performance is robust to imputation methods.

##  6 Conclusion

In the age of machine learning, the standard practice of dropping missing
cross-sectional predictors is simply untenable. Previous ML studies apply ad-
hoc imputations or data adjustments, with little discussion of their
motivation. The goal of our paper is to provide guidance on this problem.

We examine several of imputation methods in the context of several return
forecasting techniques. Our main finding is invariance: for almost all return
forecasts, all imputations lead to similar inferences about mean returns and
Sharpe ratios. This invariance happens because cross-sectional predictors are
largely independent. As a result, observed predictors are largely
uninformative about missing predictors, and all imputations lead to similar
results, including simply imputing with cross-sectional means.

More broadly, our paper adds to the evidence that the cross-section of
expected stock returns is surprisingly high-dimensional (Green, Hand and Zhang
(2014, 2017); Chen and Zimmermann (2022); Bessembinder, Burt and Hrdlicka
(2021)). An outstanding issue is the problem of missing values, which could
potentially bias these papers toward finding higher dimensionality than is
really there. We find that this literature is robust to this concern and
provide further support to the idea that one should examine dozens of
predictors to fully capture the cross-section of expected returns.

##  Appendix A Appendix

###  A.1 Extended Box-Cox Transformation

For each month-predictor combination, we symmetrically winsorize 1% of
observations and then transform the winsorized data using Hawkins and Weisberg
(2017)’s extension of the Box and Cox (1964) method. As we apply this
transformation within each month-predictor separately, this avoids look-ahead
bias in our portfolio tests (Sections 4.1 and 4.2).

Box and Cox assume the data is normally distributed after a power transform

$\displaystyle g_{BC}\left(x|\lambda\right)$
$\displaystyle=\begin{cases}\frac{x^{\lambda}-1}{\lambda}&\lambda\neq 0\\\
\log x&\lambda=0\end{cases}$

where $\lambda$ is a parameter. The common practice of taking logs to reduce
the skewness of predictors (e.g. size and B/M in Fama and French (1992)) is a
special case of this transformation ($\lambda$=0). Box and Cox recommend
estimating $\lambda$ by maximum likelihood.

A limitation of Box and Cox is that it requires that the untransformed data is
always positive. Hawkins and Weisberg (2017) apply the transformation

$\displaystyle g_{HW}\left(y|\gamma\right)$
$\displaystyle=0.5\left(y+\sqrt{y^{2}+\gamma^{2}}\right)$

and then assume $g_{BC}\left(g_{HW}\left(x|\gamma\right)|\lambda\right)$ is
normally distributed. Maximizing the likelihood then gives the estimated
parameter values $\gamma$ and $\lambda$ as well as the transformed data.

We choose the Hawkins and Weisberg transformation because it is implemented in
the popular “Companion to Applied Regression” R package (Fox and Weisberg
(2018)).

###  A.2 EM Algorithm Details

For ease of notation, we drop the $t$ subscripts.

The MVN assumption (Equation (4)) implies that the log likelihood under
complete data can be written as

$\displaystyle l\left(\mu,\Sigma|x\right)$ $\displaystyle\propto
N\log\left|\Sigma^{-1}\right|-\sum_{i=1}^{N}\text{tr}\left[\Sigma^{-1}\left(x_{i}x_{i}^{\prime}-\mu
x_{i}^{\prime}-x_{i}\mu^{\prime}-\mu\mu^{\prime}\right)\right]$

Equation (4) actually assumes $\mu=0$ but we allow for a more general $\mu$ in
this Appendix since it helps make the MVN structure more transparent. This
expression is linear in $x_{i}x_{i}^{\prime}$ and $x_{i}$, so the ignorable
quasi log-likelihood (Equation (2)) just replaces $x_{i}x_{i}^{\prime}$ and
$x_{i}$ with expected values conditional on observed data.

$\displaystyle
q\left(\mu,\Sigma|\left[x_{1}^{obs},...,x_{N}^{obs}\right]\right)$
$\displaystyle=N\log\left|\Sigma^{-1}\right|-\sum_{i=1}^{N}\text{tr}\left[\Sigma^{-1}\left(S_{i}-\mu
e_{i}^{\prime}-e_{i}\mu^{\prime}-\mu\mu^{\prime}\right)\right]$ (9)
$\displaystyle S_{i}$ $\displaystyle\equiv
E\left(X_{i}X_{i}^{\prime}|X_{i}^{obs}=x_{i}^{obs},\mu,\Sigma\right)$
$\displaystyle e_{i}$ $\displaystyle\equiv
E\left(X_{i}|X_{i}^{obs}=x_{i}^{obs},\mu,\Sigma\right)$

where $S_{i}$ and $e_{i}$ can be computed using standard conditional MVN
formulas.

To maximize Equation (9), the EM algorithm implies two steps:

  1. 1.

E-step: Given parameter guess $\mu_{old},\Sigma_{old}$, replace each
observation $i$ with the conditional expectations

$\displaystyle\hat{S}_{i}$ $\displaystyle\equiv
E\left(X_{i}X_{i}^{\prime}|X_{i}^{obs}=x_{i}^{obs},\mu_{old},\Sigma_{old}\right)$
$\displaystyle\hat{e}_{i}$ $\displaystyle\equiv
E\left(X_{i}|X_{i}^{obs}=x_{i}^{obs},\mu_{old},\Sigma_{old}\right)$

  2. 2.

M-step: Plug $\hat{S}_{i}$ and $\hat{e}_{i}$ into Equation (9) and optimize:

$\displaystyle\left(\mu_{new},\Sigma_{new}\right)$
$\displaystyle=\arg\max_{\mu,\Sigma}N\log\left|\Sigma^{-1}\right|-\sum_{i=1}^{N}\text{tr}\left[\Sigma^{-1}\left(\hat{S}_{i}-\mu\hat{e}_{i}^{\prime}-\hat{e}_{i}\mu^{\prime}-\mu\mu^{\prime}\right)\right]$

The first order conditions imply

$\displaystyle\mu_{new}$ $\displaystyle=N^{-1}\sum_{i=1}^{N}\hat{e}_{i}$
$\displaystyle\Sigma_{new}$
$\displaystyle=N^{-1}\sum_{i=1}^{N}\hat{S}_{i}-\mu_{new}\mu_{new}^{\prime}$
(10)

The algorithm ends when the maximum norm
$\max\left\\{\left\|\mu_{new}-\mu_{old}\right\|_{\infty},\left\|\Sigma_{new}-\Sigma_{old}\right\|_{\infty}\right\\}\leq
10^{-4}$.

Dempster, Laird and Rubin (1977) show that
$q\left(\mu_{new},\Sigma_{new}|\left[x_{1}^{obs},...,x_{N}^{obs}\right]\right)\geq
q\left(\mu_{old},\Sigma_{old}|\left[x_{1}^{obs},...,x_{N}^{obs}\right]\right)$,
so every iteration improves the likelihood.

###  A.3 Correlations and Eigenvalues in a Simulation

A formal relationship between the distribution of correlations and the
distribution of eigenvalues is difficult to derive. As an informal
demonstration, we simulate 125 $\times$ 125 correlation matrices using
Lewandowski, Kurowicka and Joe (2009)’s method, which models correlation
matrices by assuming partial correlations are given by a symmetric beta
distribution. This method boils down the $N$ choose $2$ parameters in a
general correlation matrix to a single parameter. The resulting distribution
of correlations is symmetric around zero, consistent with the average zero
correlations found in a variety of studies (McLean and Pontiff (2016); Green,
Hand and Zhang (2017); Chen and Zimmermann (2022)). We then conduct a PCA
variance decomposition and evaluate the imputation expression (Equation (7))
assuming that a random 1/3 of predictor pairs are missing. For simplicity, we
assume that the actual covariance matrix is known. Results average over 20
simulations.

Figure A.1 shows the result. Panel (a) shows simulations in which the beta
distribution shape parameters are both $1.2$. This setting results in a
distribution of correlations that is quite dispersed, with a sizeable share of
correlations exceeding 0.5 in absolute value (left panel). Accordingly, just 5
principal components (PCs) capture 80% of total variance (middle panel),
implying that predictor predictors are low dimensional.

[Figure A.1 about here]

Consistent with intuition, low dimensionality implies that a careful
imputation is quite different than just imputing with observed means. This
difference is seen in the right panel, which plots the distribution of
imputation slopes. More precisely, we stack all the imputation slope vectors
$\left\\{\hat{\beta}_{i,j,t}\right\\}_{i=1,...,N_{t};j=1,...,J}$ into one
large vector and plot the distribution. This distribution shows that many
slopes are meaningfully different than zero. Indeed, the mean absolute slope
is 1.7, indicating that EM imputation would be very different than simple mean
imputation.

In contrast, Panel (c) describes a high-dimensional setting. Here, the partial
correlations are given by a beta distribution with shape parameters both equal
to 15. This setting implies that almost all full correlations are less than
0.5 in absolute value, and thus even 10 PCs span less than 50% of total
variance. Consistent with intuition, the imputation slopes cluster close to
zero, with an average absolute value of 0.14. In this case, EM imputation
would be very similar to simple mean imputation.

Panel (b) shows an intermediate setting. Here, the beta shape parameters are
equal to 4, and the resulting charts show dimensionality and imputation slopes
in between those of Panels (a) and (c). Thus, the importance of careful
imputation is an issue of magnitudes, with higher dimensionality implying less
importance, in a generally smooth manner.

###  A.4 A Closer Look at EM vs Observed Correlations

Though the overall distributions are nearly identical, the EM algorithm does
have an effect on correlations at the predictor-pair level. These effects are
seen in Figure A.2, which plots the distribution of difference between EM
correlations and observed correlations for each predictor pair. In other
words, while Figure LABEL:fig:cor plots the distribution of correlations and
then compares across EM and observed, Figure A.2 compares across EM and
observed, and then plots the distribution.

[Figure A.2 about here]

Panel (a) of Figure A.2 shows that the EM and observed correlations differ on
average by about 0.05 in absolute value. This difference is tiny in absolute
terms, but given that correlations are typically close to zero the difference
is non-trivial, relatively speaking. This alternative perspective is seen in
Panel (b), which shows the substantial changes in correlations scaled by the
absolute value of the observed correlation. In previous drafts, we also showed
that EM correlations are more informative than observed correlations in the
sense that they lead to higher expected returns in PCR tests.

###  A.5 Additional Imputation Method Details

####  A.5.1 MVN EM with Time-Series Information

Our main EM imputation handles each cross-section independently. This
imputation incorporates time-series information.

We impute in three steps for each “forecasting month” $\tau$:

  1. 1.

Using OLS estimate $\hat{\phi}_{j}$ for the AR1 model

$\displaystyle X_{i,j,t}$
$\displaystyle=\phi_{j}X_{i,j,t-h}+\varepsilon_{i,j,t}$ (11)

where $X_{i,j,t}$ is Box-Cox transformed predictor (as in Section 3.2),
$h\in\left\\{1,3,12\right\\}$ is the periodicity of the underlying data
updates (based on Chen and Zimmermann (2022)’s documentation), $\phi_{j}$ is
the persistence of predictor $j$, and $\epsilon_{i,j,t}$ is a zero mean shock.
In this step, we use observed data from month $\tau-59$ to $\tau$.

  2. 2.

Using EM, estimate $\hat{\Omega}_{\tau}$ for

$\displaystyle\hat{\varepsilon}_{i,\tau}$
$\displaystyle\equiv\left[\begin{array}[]{c}\hat{\varepsilon}_{i,1,\tau}\\\
\hat{\varepsilon}_{i,2,\tau}\\\ ...\\\
\hat{\varepsilon}_{i,J,\tau}\end{array}\right]\sim\text{MVN}\left(0,\Omega_{\tau}\right).$
(16)

where $\hat{\varepsilon}_{i,j,\tau}\equiv
X_{i,j,\tau}-\hat{\phi}_{j}X_{i,j,t-h}$. Assuming standard regularity
conditions, this two-step estimation is consistent (Wooldridge (1994)).

  3. 3.

For missing $X_{i,j,\tau}$, impute by combining the prediction of the AR1
model with the residual estimate—that is, use
$\hat{\phi}_{j}X_{i,j,t-h}+\hat{\varepsilon}_{i,j,\tau}$ (if $X_{i,j,t-h}$ is
observed) or $\hat{\varepsilon}_{i,j,\tau}$ (if $X_{i,j,t-h}$ is not
observed).

####  A.5.2 Practical EM for Probabilistic PCA

We did not restrict the covariance matrix in our main EM imputation. This
imputation imposes a factor structure and estimates using a variant of EM
(Roweis (1997), see also Tipping and Bishop (1999) ).

We estimate the following model:

$\displaystyle X_{i,t}$ $\displaystyle=\Lambda_{t}F_{i,t}+\varepsilon_{i,t}$
(17) $\displaystyle F_{i,t}$ $\displaystyle\sim\text{MVN}\left(0,I\right)$
$\displaystyle\varepsilon_{i,t}$
$\displaystyle\sim\text{MVN}\left(0,\sigma^{2}I\right)$ (18)

where $X_{i,t}$ is a vector of predictors for stock $i$ in month $t$,
$F_{i,t}$ is a $K$-dimensional vector of factor realizations for stock $i$ in
month $t$, $\Lambda_{t}$ is a $N_{t}\times K$ matrix that maps factors into
predictors, and $\sigma^{2}$ is a scalar.

The estimation alternates between two steps for each forecasting month $t$:

  1. 1.

E-step: Given a guess for $\Lambda_{t}$ and $F_{i,t}$, replace the missing
components of $X_{i,t}$ with $\Lambda_{t}F_{i,t}$. Find the expectations of
$F_{i,t}$ and $F_{i,t}F_{i,t}^{\prime}$ conditional on $\Lambda_{t}$,
$\sigma$, and $X_{i,t}$.

  2. 2.

M-step: Plug the expectations of $F_{i,t}$ and $F_{i,t}F_{i,t}^{\prime}$ into
the log-likelihood assuming that all $X_{i,t}$ and $F_{i,t}$ are observed, and
maximize to obtain new $\Lambda_{t}$ and $F_{i,t}$.

The algorithm iterates these two steps until the objective in the M-step no
longer improves.

As in the case of our baseline model (Equation (4)), Equations (17)-(18) imply
that the ignorable quasi-likelihood (Equation (2)) can be optimized using EM
(Tipping and Bishop (1999)). This algorithm does not follow the full EM
method, as it does not include higher moments of missing $X_{i,t}$ in the
E-step. However, the computer science literature finds that less formal
variants of EM work well on this model using simulated data (Roweis (1997);
Ilin and Raiko (2010)). Our implementation uses the R library pcaMethods,
which follows Roweis (1997)’s algorithm.

###  A.6 Additional Forecasting Method Details

GBRT details are fully described in Sections 5.2 and 5.3.

Neural network details follow Gu et al. (2020). NN1 has a single hidden layer
of 32 neurons and NN3 has hidden layers with 32, 16, and 8 neurons. All layers
are connected with the ReLU activiation function. The objective is L2 with a
L1 penalty to weight parameters and is minimized using the Adam extension of
Stochastic Gradient Descent under early stopping with a patience parameter of
5 in batches of 10,000, for 100 epochs. The L1 penalty and learning rate are
tuned as described in the main text. Last, we form final return forecasts as
the ensemble average of 10 neural network forecasts.

###  A.7 Appendix Figures

Figure A.1: Dimensionality and Imputation Simulations

We simulate correlation matricies, compute eigenvalue decompositions, and
imputation slopes. Imputation slopes are the response of an imputed missing
value to a one standard deviation change in an observed predictor (Equation
(7)). Slopes are the elements of the vectors $\hat{\beta}_{i,j,t}^{\prime}$
computed with Equation (8), reprinted here:

$\displaystyle\hat{\beta}_{i,j,t}^{\prime}$
$\displaystyle=\hat{\Sigma}_{j,obs|i,t}\hat{\Sigma}_{obs|i,obs|i,t}^{-1}.$

Simulations assume symmetric-beta-distributed partial correlations following
Lewandowski et al. (2009). Panel (a) uses a beta parameter of 1.2, Panel (b)
uses 4, and Panel (c) uses 15. Interpretation: High dimensionality implies
that slopes are close to zero, which implies simple mean imputation is similar
to EM imputation.

(a)

(b)

(c)

Figure A.2: A Closer Look at EM vs Observed Correlations

Data is from June 2000. We compute pairwise correlations using the EM
algorithm and observed available cases, and then take differences (Panel (a)).
Panel (b) normalizes the differences by the absolute value of the observed
available case correlation. Interpretation: Though the distribution of
correlations is visually identical using EM or observed available case (Figure
LABEL:fig:cor), at the predictor-pair level the correlations differ, and with
the median difference about 30% of the available case correlation.

(a)

(b)

Figure A.3: BLLP’s Imputation and Scaled-PCA Returns

Each month, we forecast returns using scaled principal components regression
and form long-short deciles portfolios. ‘BLLP loc B-XS’ imputes using the
local backward cross-sectional model in BLLP, where we use six factors to be
consistent with Section 4.1 of BLLP, and lags used in the second stage
regression are restricted to relevant predictor update months based on
documentation for the Chen and Zimmermann (2022) dataset. Interpretation:
Taking on time-series information following Bryzgalova et al. (2023) results
in similar estimates as simple mean imputation.

(a)

(b)

## References

  * Azen et al. (1989) Azen, S.P., van Guilder, M., Hill, M.A., 1989.  Estimation of parameters and missing values under a regression model with non-normally distributed and non-randomly incomplete data.  Statistics in Medicine 8, 217–228.
  * Azevedo et al. (2022) Azevedo, V., Kaiser, S., Müller, S., 2022.  Stock market anomalies and machine learning across the globe.  Available at SSRN 4071852 .
  * Beckmeyer and Wiedemann (2022) Beckmeyer, H., Wiedemann, T., 2022\.  Recovering missing firm characteristics with attention-based machine learning.  Available at SSRN 4003455 .
  * Bessembinder et al. (2021) Bessembinder, H., Burt, A., Hrdlicka, C.M., 2021.  Time series variation in the factor zoo.  Aaron Paul and Hrdlicka, Christopher M., Time Series Variation in the Factor Zoo (December 22, 2021) .
  * Bowles et al. (2020) Bowles, B., Reed, A.V., Ringgenberg, M.C., Thornock, J.R., 2020\.  Anomaly time.  Available at SSRN 3069026 .
  * Box and Cox (1964) Box, G.E., Cox, D.R., 1964.  An analysis of transformations.  Journal of the Royal Statistical Society: Series B (Methodological) 26, 211–243.
  * Bryzgalova et al. (2023) Bryzgalova, S., Lerner, S., Lettau, M., Pelger, M., 2023\.  Missing financial data.  Available at SSRN 4106794 .
  * Chen and Velikov (2022) Chen, A.Y., Velikov, M., 2022\.  Zeroing in on the expected returns of anomalies.  Journal of Financial and Quantitative Analysis , 1–83.
  * Chen and Zimmermann (2022) Chen, A.Y., Zimmermann, T., 2022\.  Open source cross sectional asset pricing.  Critical Finance Review .
  * Chen et al. (2023) Chen, L., Pelger, M., Zhu, J., 2023.  Deep learning in asset pricing.  Management Science .
  * Da et al. (2022) Da, R., Nagel, S., Xiu, D., 2022.  The Statistical Limit of Arbitrage.  Technical Report. Technical Report, Chicago Booth.
  * DeMiguel et al. (2020) DeMiguel, V., Martin-Utrera, A., Nogales, F.J., Uppal, R., 2020\.  A transaction-cost perspective on the multitude of firm characteristics.  The Review of Financial Studies 33, 2180–2222.
  * Dempster et al. (1977) Dempster, A.P., Laird, N.M., Rubin, D.B., 1977.  Maximum likelihood from incomplete data via the em algorithm.  Journal of the Royal Statistical Society: Series B (Methodological) 39, 1–22.
  * Efron and Hastie (2016) Efron, B., Hastie, T., 2016\.  Computer age statistical inference. volume 5.  Cambridge University Press.
  * Engelberg et al. (2018) Engelberg, J., McLean, R.D., Pontiff, J., 2018.  Anomalies and news.  The Journal of Finance 73, 1971–2001.
  * Fama and French (1992) Fama, E., French, K., 1992.  The cross-section of expected stock returns.  Journal of finance , 427–465.
  * Fama and French (1993) Fama, E., French, K., 1993.  Common risk factors in the returns on stocks and bonds.  Journal of financial economics 33, 3–56.
  * Fama and French (2015) Fama, E.F., French, K.R., 2015\.  A five-factor asset pricing model.  Journal of Financial Economics 116, 1 – 22.  URL: http://www.sciencedirect.com/science/article/pii/S0304405X14002323, doi:https://doi.org/10.1016/j.jfineco.2014.10.010.
  * Fox and Weisberg (2018) Fox, J., Weisberg, S., 2018\.  An R companion to applied regression.  Sage publications.
  * Freyberger et al. (2022) Freyberger, J., Höppner, B., Neuhierl, A., Weber, M., 2022\.  Missing data in asset pricing panels.  University of Chicago Becker Friedman Institute Working Paper .
  * Freyberger et al. (2020) Freyberger, J., Neuhierl, A., Weber, M., 2020.  Dissecting characteristics nonparametrically.  The Review of Financial Studies 33, 2326–2377.
  * Fulop et al. (2022) Fulop, A., Li, J., Wang, M., 2022.  Option mispricing and alpha portfolios.  Available at SSRN .
  * Graham (2009) Graham, J.W., 2009.  Missing data analysis: Making it work in the real world.  Annual review of psychology 60, 549–576.
  * Graham and Schafer (1999) Graham, J.W., Schafer, J.L., 1999\.  On the performance of multiple imputation for multivariate data with small sample size.  Statistical strategies for small sample research 50, 1–27.
  * Green et al. (2014) Green, J., Hand, J., Zhang, F., 2014.  The remarkable multidimensionality in the cross-section of expected us stock returns.  Available at SSRN 2262374.
  * Green et al. (2013) Green, J., Hand, J.R., Zhang, X.F., 2013.  The supraview of return predictive signals.  Review of Accounting Studies 18, 692–730.
  * Green et al. (2017) Green, J., Hand, J.R., Zhang, X.F., 2017.  The characteristics that provide independent information about average us monthly stock returns.  The Review of Financial Studies , hhx019.
  * Gu et al. (2020) Gu, S., Kelly, B., Xiu, D., 2020.  Empirical asset pricing via machine learning.  The Review of Financial Studies 33, 2223–2273.
  * Haitovsky (1968) Haitovsky, Y., 1968.  Missing data in regression analysis.  Journal of the Royal Statistical Society: Series B (Methodological) 30, 67–82.
  * Han et al. (2021) Han, Y., He, A., Rapach, D., Zhou, G., 2021.  Expected stock returns and firm characteristics: E-lasso, assessment, and implications.  Assessment, and Implications (September 10, 2021) .
  * Han et al. (2022) Han, Y., He, A., Rapach, D., Zhou, G., 2022.  Expected stock returns and firm characteristics: E-enet, assessment, and implications.  Working Paper .
  * Hartley (1958) Hartley, H.O., 1958.  Maximum likelihood estimation from incomplete data.  Biometrics 14, 174–194.
  * Harvey et al. (2016) Harvey, C.R., Liu, Y., Zhu, H., 2016.  … and the cross-section of expected returns.  The Review of Financial Studies 29, 5–68.
  * Haugen and Baker (1996) Haugen, R.A., Baker, N.L., 1996\.  Commonality in the determinants of expected stock returns.  Journal of Financial Economics 41, 401–439.
  * Hawkins and Weisberg (2017) Hawkins, D., Weisberg, S., 2017\.  Combining the box-cox power and generalised log transformations to accommodate nonpositive responses in linear and mixed-effects linear models.  South African Statistical Journal 51, 317–328.
  * Huang et al. (2022) Huang, D., Jiang, F., Li, K., Tong, G., Zhou, G., 2022\.  Scaled pca: A new approach to dimension reduction.  Management Science 68, 1678–1695.
  * Ilin and Raiko (2010) Ilin, A., Raiko, T., 2010.  Practical approaches to principal component analysis in the presence of missing values.  The Journal of Machine Learning Research 11, 1957–2000.
  * Jensen et al. (2022) Jensen, T.I., Kelly, B.T., Malamud, S., Pedersen, L.H., 2022\.  Machine learning and the implementable efficient frontier.  Available at SSRN 4187217 .
  * Jensen et al. (Forthcoming) Jensen, T.I., Kelly, B.T., Pedersen, L.H., Forthcoming.  Is there a replication crisis in finance?  Journal of Finance .
  * Jolliffe (2002) Jolliffe, I.T., 2002.  Principal Component Analysis.  Springer.
  * Kelly et al. (2023) Kelly, B., Malamud, S., Pedersen, L.H., 2023.  Principal portfolios.  The Journal of Finance .
  * Kozak et al. (2020) Kozak, S., Nagel, S., Santosh, S., 2020.  Shrinking the cross-section.  Journal of Financial Economics 135, 271–292.
  * Lettau and Pelger (2020) Lettau, M., Pelger, M., 2020\.  Factors that fit the time series and cross-section of stock returns.  The Review of Financial Studies 33, 2274–2325.
  * Lewandowski et al. (2009) Lewandowski, D., Kurowicka, D., Joe, H., 2009.  Generating random correlation matrices based on vines and extended onion method.  Journal of multivariate analysis 100, 1989–2001.
  * Lewellen (2015) Lewellen, J., 2015.  The cross-section of expected stock returns.  Critical Finance Review 4, 1–44.
  * Li et al. (2022) Li, B., Rossi, A., Yan, S., Zheng, L., 2022.  Real-time machine learning in the cross-section of stock returns.  Working Paper.
  * Little (1988) Little, R.J., 1988.  Robust estimation of the mean and covariance matrix from data with missing values.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 37, 23–38.
  * Little (1992) Little, R.J., 1992.  Regression with missing x’s: a review.  Journal of the American statistical association 87, 1227–1237.
  * Little (2021) Little, R.J., 2021.  Missing data assumptions.  Annual Review of Statistics and Its Application 8, 89–107.
  * Little and Rubin (2019) Little, R.J., Rubin, D.B., 2019\.  Statistical analysis with missing data. volume 793\.  Wiley.
  * Lopez-Lira and Roussanov (2020) Lopez-Lira, A., Roussanov, N.L., 2020\.  Do common factors really explain the cross-section of stock returns?  Available at SSRN 3628120 .
  * McLean and Pontiff (2016) McLean, R.D., Pontiff, J., 2016\.  Does academic research destroy stock return predictability?  The Journal of Finance 71, 5–32.
  * Orchard and Woodbury (1972) Orchard, T., Woodbury, M.A., 1972\.  A missing information principle: theory and applications, in: Volume 1 Theory of Statistics. University of California Press, pp. 697–716.
  * Roweis (1997) Roweis, S., 1997.  Em algorithms for pca and spca.  Advances in neural information processing systems 10\.
  * Rubin (1976) Rubin, D.B., 1976.  Inference and missing data.  Biometrika 63, 581–592.
  * Simon et al. (2022) Simon, F., Weibels, S., Zimmermann, T., 2022.  Deep parametric portfolio policies.  Available at SSRN 4150292 .
  * Tipping and Bishop (1999) Tipping, M.E., Bishop, C.M., 1999\.  Probabilistic principal component analysis.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 61, 611–622.
  * Titman et al. (2004) Titman, S., Wei, K.J., Xie, F., 2004.  Capital investments and stock returns.  Journal of financial and Quantitative Analysis 39, 677–700.
  * Wooldridge (1994) Wooldridge, J.M., 1994.  Estimation and inference for dependent processes.  Handbook of econometrics 4, 2639–2738.
  * Xiong and Pelger (Forthcoming) Xiong, R., Pelger, M., Forthcoming.  Large dimensional latent factor modeling with missing observations and applications to causal inference.  Journal of Econometrics .
  * Zhang and Aaraba (2022) Zhang, C., Aaraba, A., 2022\.  Dynamic optimal portfolio with trading cost constraints via deep reinforcement learning.  Available at SSRN .

Table 1:  Predictors Ranked by Observations in 1985: Data are 125 predictors from Chen and Zimmermann (2022) with the most observations in 1985, are continuous, and have at least 2 observations in every month between 1985 and 2020. ‘% Observed’ is the average share of stocks that have data across months in 1985. Interpretation: This table provides a guide to the origins of missingness that are discussed in Section 2.2.  Acronym | Reference | % Observed | Rank
---|---|---|---
streversal | Jegadeesh 1989 | 100.0 | 1
maxret | Bali, Cakici, and Whitelaw 2010 | 98.3 | 2
idiorisk | Ang et al. 2006 | 97.7 | 3
idiovol3f | Ang et al. 2006 | 97.6 | 4
returnskew3f | Bali, Engle and Murray 2015 | 97.6 | 5
price | Blume and Husic 1972 | 97.5 | 6
size | Banz 1981 | 97.5 | 7
high52 | George and Hwang 2004 | 97.5 | 8
zerotradealt1 | Liu 2006 | 97.2 | 9
dolvol | Brennan, Chordia, Subra 1998 | 96.2 | 10
idiovolaht | Ali, Hwang, and Trombley 2003 | 95.0 | 11
returnskew | Bali, Engle and Murray 2015 | 93.8 | 12
zerotrade | Liu 2006 | 93.3 | 13
indmom | Grinblatt and Moskowitz 1999 | 92.7 | 14
coskewness | Harvey and Siddique 2000 | 90.9 | 15
volmkt | Haugen and Baker 1996 | 90.7 | 16
pricedelayrsq | Hou and Moskowitz 2005 | 90.6 | 17
pricedelayslope | Hou and Moskowitz 2005 | 90.6 | 18
mom6m | Jegadeesh and Titman 1993 | 89.4 | 19
mom12moffseason | Heston and Sadka 2008 | 88.7 | 20
illiquidity | Amihud 2002 | 88.6 | 21
zerotradealt12 | Liu 2006 | 88.1 | 22
pricedelaytstat | Hou and Moskowitz 2005 | 87.2 | 23
momseasonshort | Heston and Sadka 2008 | 85.6 | 24
mom12m | Jegadeesh and Titman 1993 | 85.3 | 25
intmom | Novy-Marx 2012 | 84.8 | 26
beta | Fama and MacBeth 1973 | 81.4 | 27
mrreversal | De Bondt and Thaler 1985 | 79.6 | 28
shareiss1y | Pontiff and Woodgate 2008 | 79.5 | 29
herf | Hou and Robinson 2006 | 79.1 | 30
volsd | Chordia, Subra, Anshuman 2001 | 77.0 | 31
bookleverage | Fama and French 1992 | 76.2 | 32
roe | Haugen and Baker 1996 | 76.2 | 33
opleverage | Novy-Marx 2010 | 75.9 | 34
hire | Bazdresch, Belo and Lin 2014 | 75.6 | 35
momseason | Heston and Sadka 2008 | 72.8 | 36
momoffseason | Heston and Sadka 2008 | 72.7 | 37
am | Fama and French 1992 | 72.4 | 38
cf | Lakonishok, Shleifer, Vishny 1994 | 72.4 | 39
sp | Barbee, Mukherji and Raines 1996 | 72.1 | 40
leverage | Bhandari 1988 | 72.0 | 41
cashprod | Chandrashekar and Rao 2009 | 71.8 | 42
bpebm | Penman, Richardson and Tuna 2007 | 71.7 | 43
ebm | Penman, Richardson and Tuna 2007 | 71.7 | 44
bm | Rosenberg, Reid, and Lanstein 1985 | 71.2 | 45
numearnincrease | Loh and Warachka 2012 | 71.2 | 46
tax | Lev and Nissim 2004 | 69.0 | 47
bidaskspread | Amihud and Mendelsohn 1986 | 67.8 | 48
xfin | Bradshaw, Richardson, Sloan 2006 | 67.3 | 49
assetgrowth | Cooper, Gulen and Schill 2008 | 66.9 | 50
chinv | Thomas and Zhang 2002 | 66.9 | 51
delcoa | Richardson et al. 2005 | 66.9 | 52
dellti | Richardson et al. 2005 | 66.9 | 53
delequ | Richardson et al. 2005 | 66.8 | 54
dnoa | Hirshleifer, Hou, Teoh, Zhang 2004 | 66.8 | 55
betafp | Frazzini and Pedersen 2014 | 66.7 | 56
volumetrend | Haugen and Baker 1996 | 66.7 | 57
gp | Novy-Marx 2013 | 66.6 | 58
grltnoa | Fairfield, Whisenant and Yohn 2003 | 66.6 | 59
netequityfinance | Bradshaw, Richardson, Sloan 2006 | 66.6 | 60
totalaccruals | Richardson et al. 2005 | 66.5 | 61
noa | Hirshleifer et al. 2004 | 66.5 | 62
equityduration | Dechow, Sloan and Soliman 2004 | 66.4 | 63
accruals | Sloan 1996 | 66.4 | 64
chnwc | Soliman 2008 | 66.4 | 65
delcol | Richardson et al. 2005 | 66.4 | 66
chnncoa | Soliman 2008 | 66.2 | 67
delfinl | Richardson et al. 2005 | 66.2 | 68
delnetfin | Richardson et al. 2005 | 66.2 | 69
bmdec | Fama and French 1992 | 66.2 | 70
cheq | Lockwood and Prombutr 2010 | 65.0 | 71
earnsupbig | Hou 2007 | 64.4 | 72
indretbig | Hou 2007 | 64.4 | 73
chtax | Thomas and Zhang 2011 | 63.6 | 74
rds | Landsman et al. 2011 | 63.4 | 75
roaq | Balakrishnan, Bartov and Faurel 2010 | 63.1 | 76
investppeinv | Lyandres, Sun and Zhang 2008 | 62.7 | 77
chinvia | Abarbanell and Bushee 1998 | 62.4 | 78
pctacc | Hafzalla, Lundholm, Van Winkle 2011 | 62.0 | 79
residualmomentum | Blitz, Huij and Martens 2011 | 61.9 | 80
lrreversal | De Bondt and Thaler 1985 | 61.9 | 81
abnormalaccruals | Xie 2001 | 61.8 | 82
announcementreturn | Chan, Jegadeesh and Lakonishok 1996 | 61.2 | 83
cboperprof | Ball et al. 2016 | 61.0 | 84
herfasset | Hou and Robinson 2006 | 60.5 | 85
herfbe | Hou and Robinson 2006 | 60.5 | 86
netdebtfinance | Bradshaw, Richardson, Sloan 2006 | 58.8 | 87
cfp | Desai, Rajgopal, Venkatachalam 2004 | 58.6 | 88
entmult | Loughran and Wellman 2011 | 58.2 | 89
varcf | Haugen and Baker 1996 | 56.9 | 90
grsaletogrinv | Abarbanell and Bushee 1998 | 55.8 | 91
grsaletogroverhead | Abarbanell and Bushee 1998 | 55.2 | 92
operprofrd | Ball et al. 2016 | 55.1 | 93
grcapx | Anderson and Garcia-Feijoo 2006 | 54.9 | 94
ep | Basu 1977 | 54.5 | 95
earningssurprise | Foster, Olsen and Shevlin 1984 | 54.2 | 96
chassetturnover | Soliman 2008 | 52.7 | 97
revenuesurprise | Jegadeesh and Livnat 2006 | 51.5 | 98
shareiss5y | Daniel and Titman 2006 | 51.4 | 99
cash | Palazzo 2012 | 50.3 | 100
grcapx3y | Anderson and Garcia-Feijoo 2006 | 50.1 | 101
momseason06yrplus | Heston and Sadka 2008 | 49.0 | 102
momoffseason06yrplus | Heston and Sadka 2008 | 48.7 | 103
betatailrisk | Kelly and Jiang 2014 | 47.6 | 104
invgrowth | Belo and Lin 2012 | 45.9 | 105
investment | Titman, Wei and Xie 2004 | 43.9 | 106
feps | Cen, Wei, and Zhang 2006 | 43.7 | 107
analystrevision | Hawkins, Chamberlin, Daniel 1984 | 43.1 | 108
momseason11yrplus | Heston and Sadka 2008 | 42.0 | 109
momoffseason11yrplus | Heston and Sadka 2008 | 41.9 | 110
std_turn | Chordia, Subra, Anshuman 2001 | 41.7 | 111
intancfp | Daniel and Titman 2006 | 41.2 | 112
intanep | Daniel and Titman 2006 | 41.2 | 113
compositedebtissuance | Lyandres, Sun and Zhang 2008 | 41.1 | 114
intansp | Daniel and Titman 2006 | 41.1 | 115
netpayoutyield | Boudoukh et al. 2007 | 40.7 | 116
intanbm | Daniel and Titman 2006 | 40.4 | 117
exclexp | Doyle, Lundholm and Soliman 2003 | 39.6 | 118
meanrankrevgrowth | Lakonishok, Shleifer, Vishny 1994 | 39.5 | 119
netdebtprice | Penman, Richardson and Tuna 2007 | 38.3 | 120
compequiss | Daniel and Titman 2006 | 38.1 | 121
trendfactor | Han, Zhou, Zhu 2016 | 37.1 | 122
operprof | Fama and French 2006 | 35.9 | 123
tang | Hahn and Lee 2009 | 34.0 | 124
payoutyield | Boudoukh et al. 2007 | 33.4 |
Table 1: (continued)

Data are 125 predictors from Chen and Zimmermann (2022) with the most
observations in 1985, are continuous, and have at least 2 observations in
every month between 1985 and 2020. Panel (a) computes the % of stocks with
observations for individual predictors and then computes percentiles across
predictors. Panel (b) selects the $J$ most-observed predictors (Table
LABEL:tab:final-list) and then computes the % observed among these using two
methods: ‘Predictor-Stocks’ sums predictor-stock observations and then divides
by $J\times$ the number of stocks. ‘Stocks w/ all $J$’ counts the number of
stocks with all $J$ predictors and divides by the number of stocks.
Interpretation: Missing values are a much larger problem when combining many
predictors. When combining 50 or more predictors, one must impute missing
values or drop most stocks.

Panel (a): Individual Predictors
---
Predictor |  | % of Stocks With Predictor Data in June |
Percentile |  | 1970 | 1975 | 1980 | 1985 | 1990 | 1995 | 2000
5 |  | 0.0 | 24.8 | 37.5 | 39.6 | 36.1 | 36.8 | 39.2
25 |  | 56.8 | 43.6 | 54.4 | 54.9 | 54.3 | 53.2 | 56.1
50 |  | 69.8 | 63.2 | 69.1 | 66.4 | 68.2 | 70.8 | 71.8
75 |  | 80.6 | 68.7 | 77.7 | 76.2 | 77.7 | 76.9 | 75.8
95 |  | 98.3 | 94.0 | 96.0 | 97.5 | 97.8 | 97.7 | 97.1

Panel (b): Combining $J$ Predictors
---
|  | June 1975 |  | June 1985 |  | June 1995
|  | Predictor- |  Stocks w/ $J$ |  | Predictor- |  Stocks w/ $J$ |  | Predictor- |  Stocks w/ $J$
$J$ |  | Stocks w/ | Predictors |  | Stocks w/ | Predictors |  | Stocks w/ | Predictors
|  | Data (%) | (%) |  | Data (%) | (%) |  | Data (%) | (%)
25 |  | 83.1 | 45.8 |  | 93.1 | 77.5 |  | 91.6 | 72.0
50 |  | 77.5 | 32.9 |  | 83.5 | 35.2 |  | 83.6 | 44.6
75 |  | 71.9 | 0.0 |  | 77.7 | 15.7 |  | 78.8 | 21.7
100 |  | 66.9 | 0.0 |  | 72.7 | 7.6 |  | 73.7 | 8.3
125 |  | 59.5 | 0.0 |  | 66.5 | 0.2 |  | 66.9 | 0.2

Strategies long the 500 stocks with the most positive predictors and short 500
stocks with the most negative predictors. ‘EM Imputed’ uses multivariate
normal imputation (Section 3.2) and ‘Missing Dropped’ drops stock-months with
missing predictors. For both methods, we drop predictor-months with less than
1000 stock-predictor observations. _‘EM Improvement’_ is ‘EM Imputed’ minus
‘Missing Dropped,’ computed at the predictor level. Interpretation: Though
careful handling of missing values generally improves performance, the effect
overall is small.

|  | Percentile Across Predictors
---|---|---
|  | 10 | 25 | 50 | 75 | 90
Panel (a): Equal-Weighted
Mean Return (% Ann)
| EM Imputed | 0.1 | 3.4 | 7.0 | 11.5 | 15.8
| Missing Dropped | -0.3 | 2.8 | 5.9 | 10.0 | 15.1
| EM Improvement | -1.0 | -0.1 | 0.2 | 1.3 | 3.6
Sharpe Ratio (% Ann)
| EM Imputed | 0.02 | 0.23 | 0.56 | 0.90 | 1.27
| Missing Dropped | -0.02 | 0.22 | 0.53 | 0.86 | 1.13
| EM Improvement | -0.08 | -0.03 | 0.00 | 0.04 | 0.17
Panel (b): Value-Weighted
Mean Return (% Ann)
| EM Imputed | -3.0 | 0.7 | 2.7 | 6.6 | 10.9
| Missing Dropped | -2.7 | 0.5 | 2.4 | 5.9 | 9.3
| EM Improvement | -0.9 | -0.2 | 0.1 | 0.9 | 1.9
Sharpe Ratio (% Ann)
| EM Imputed | -0.15 | 0.04 | 0.20 | 0.43 | 0.60
| Missing Dropped | -0.14 | 0.03 | 0.19 | 0.36 | 0.53
| EM Improvement | -0.06 | -0.02 | 0.00 | 0.05 | 0.11

We impute missing values using one of six methods, forecast returns using one
of six methods, and then form long-short decile portfolios. Imputations,
forecasts and portfolios are formed recursively each June between 1995-2020.
‘Simple Mean’ imputes with cross-sectional means, ‘MVN EM’ is maximum
likelihood under the multivariate normal model (Section 3.2). Other imputation
definitions are found in Section 5.1. ‘PCR’ and ‘sPCR’ run OLS on PCs or
scaled PCs (Section 4.2). Other forecasts are defined in Section 5.2.
Hyperparameter tuning is described in Section 5.3. OLS is marked with †
because it uses no tuning. Interpretation: Inferences vary little across
imputations, including ones that use time-series information. OLS performs
well, consistent with the high-dimensionality of expected returns.

Panel (b): Value-Weighted Mean Return (% Annualized) Forecasting Methods
Linear Non-Linear OLS† PCR sPCR GBRT NN3 NN1 _Ad-hoc Imputations_ Simple Mean
21 17 20 14 35 31 Industry / Size Dec 22 21 19 17 36 30 Last Observed 20 19 20
9 36 31 _EM Imputations_ MVN EM 25 20 23 20 33 27 MVN EM on AR1 Resid 19 16 20
18 30 29 Practical EM for pPCA 23 18 20 19 33 30

Data are the most-observed 125 predictors in 1985. Shaded (red) areas indicate
missing values. Interpretation: When combining 125 predictors, a stock can
only be used without imputation if it has values for all 125 predictors,
corresponding to column that is entirely non-shaded. Almost no columns satisfy
this requirement. Missingness patterns are complex, suggesting that a method
that handles general patterns, like the EM algorithm, is required.

Each month, we forecast next month’s returns using PCR on the past 10 years of
data. Strategies go long stocks in the top decile of forecasted returns and
short the bottom decile. ‘EM’ imputes missing values with the EM algorithm.
‘Simple mean’ imputes with cross-sectional means. Predictor data runs
1985-2020. Real-time returns use data from 1995-2020. Interpretation: Simple
mean imputation works surprisingly well. Many PCs are required to capture the
potential mean returns.

Each month, we form long-short portfolios based scaled-PCA forecasts (Huang et
al. (2022)) using the past 10 years of data. ‘EM’ imputes missing values with
the EM algorithm. ‘Simple mean’ imputes with cross-sectional means. Predictor
data runs 1985-2020. Real-time returns use data from 1995-2020.
Interpretation: Scaled-PCs lead to higher mean returns than standard PCs
(Figure LABEL:fig:pcr). Nevertheless, the results still indicate many
dimensions of predictability. Simple mean and EM imputation lead to similar
inferences.

All predictors are transformed using an extended Box-Cox algorithm to have
approximately normal marginal distributions at the predictor level (Section
3.2). ‘Observed’ uses all available cases at the predictor-pair level (a.k.a.
pairwise-complete). ‘EM Algo’ uses the EM-estimated covariance matrix
($\Sigma_{t}$ from Equation (4)). Interpretation: The correlations cluster
near zero, implying that imputation slopes are close to zero (Equation (8)),
and thus simple mean imputation should work fairly well.

Imputation slopes are the response of an imputed missing value to a one
standard deviation change in an observed predictor (Equation (7)). Slopes are
the elements of the vector $\hat{\beta}_{i,j,t}^{\prime}$ from Equation (8),
reprinted here:

$\displaystyle\hat{\beta}_{i,j,t}^{\prime}$
$\displaystyle=\hat{\Sigma}_{j,obs|i,t}\hat{\Sigma}_{obs|i,obs|i,t}^{-1}.$

Interpretation: Consistent with Figure LABEL:fig:cor, the high dimensionality
of predictors implies imputation slopes are close to zero, and that simple
mean and EM imputation are similar.