# Bayesian Imputation with Optimal Look-Ahead-Bias and Variance
Tradeoff††thanks: We gratefully acknowledge the financial support and feedback
of MSCI and its data science team. Material in this paper is based upon work
supported by the Air Force Office of Scientific Research under award number
FA9550-20-1-0397. Additional support is gratefully acknowledged from NSF
grants 1915967, 1820942, 1838676, and also from the China Merchant Bank.

Jose Blanchet Stanford University, Department of Management Science &
Engineering, Email<EMAIL_ADDRESS>Fernando Hernandez Stanford
University, Department of Management Science & Engineering, Email:
<EMAIL_ADDRESS>Viet Anh Nguyen Stanford University, Department of
Management Science & Engineering, Email<EMAIL_ADDRESS>Markus
Pelger Stanford University, Department of Management Science & Engineering,
Email<EMAIL_ADDRESS>Xuhui Zhang Stanford University, Department of
Management Science & Engineering, Email<EMAIL_ADDRESS>
(January 31, 2022)

Missing time-series data is a prevalent problem in finance. Imputation methods
for time-series data are usually applied to the full panel data with the
purpose of training a model for a downstream out-of-sample task. For example,
the imputation of missing returns may be applied prior to estimating a
portfolio optimization model. However, this practice can result in a look-
ahead-bias in the future performance of the downstream task. There is an
inherent trade-off between the look-ahead-bias of using the full data set for
imputation and the larger variance in the imputation from using only the
training data. By connecting layers of information revealed in time, we
propose a Bayesian consensus posterior that fuses an arbitrary number of
posteriors to optimally control the variance and look-ahead-bias trade-off in
the imputation. We derive tractable two-step optimization procedures for
finding the optimal consensus posterior, with Kullback-Leibler divergence and
Wasserstein distance as the measure of dissimilarity between posterior
distributions. We demonstrate in simulations and an empirical study the
benefit of our imputation mechanism for portfolio optimization with missing
returns.

Keywords: Imputation, missing data, look-ahead-bias, Bayesian consensus
posterior, Kullback-Leiber divergence, Wasserstein distance, portfolio
optimization

JEL classification: C02, C11, C22, G11

##  1 Introduction

Missing time-series data is a prevalent problem in finance. Its role is
becoming even more critical due to the vast amount of data that is being
collected nowadays and the reliance on data driven analytics. Financial data
can be missing for a variety of reasons, which include for example infrequent
trading due to illiquidity, mixed frequency reporting or improper data
collection and storage. How this missing data is handled has substantial
impact on any down-stream task in risk management. A standard approach is to
impute missing data. There is a growing body of literature (Little and Rubin,
2002; van Buuren, 2012) and over 150 practical implementations available
(Mayer, Josse, Tierney, and Vialaneix, 2019) that handle the issue of imputing
missing data. However, the implications of imputation for a downstream out-of-
sample task have been neglected so far.

This paper studies the imputation of missing data in a panel which is
subsequently used for training a model for a downstream out-of-sample task.
Our leading example is the imputation of missing returns prior to estimating a
portfolio optimization model. Imputation methods are usually applied to the
full panel to impute missing values. However, this practice can result in a
look-ahead-bias in the future performance of the downstream task. Typically, a
machine learning or portfolio optimization down-stream task splits the panel
data into a train-validation-test set according to the time-ordering of
underlying events. Training corresponds to the left-most block portion of the
panel, validation corresponds to a block in the middle of the panel, and
testing is done in a block at the right end of the panel. The “globally
observed data” is the entire available panel data set representing all the
data collected during the entire time horizon. Meanwhile, “locally observed
data” refers in our discussion to data available up to a certain time in the
past (e.g. the training section of the panel data set). A standard imputation
procedure in practice is to apply a given imputation tool at hand using the
globally observed data. While being natural for the purpose of reducing error
in the imputation, ignoring the time-order of a future downstream task poses a
problem, which is of particular importance for portfolio optimization.
Portfolio selection procedures are designed to exploit subtle signals that can
lead to favorable future investment outcomes. A global imputation is likely to
create a spurious signal due to the correlation between training and testing
data induced by the imputation. The training portion of the downstream task
itself can (and often will) pickup the spurious signal created by the
imputation procedure, thus giving a false sense of good performance. For
example, if the mean return of an asset is large on the testing data, a global
imputation might increase the estimated mean on the training data and bias the
stock selection. This is an example of an overfitting phenomenon created by
“looking into the future” which investors have to mitigate. In summary, even
if the portfolio analysis is done “properly out-of-sample” on the full panel
of imputed values, the look-ahead-bias from the imputation can render the
results invalid.

Our paper provides a new conceptual perspective by connecting the imputation
with a downstream out-of-sample task. There is an inherent trade-off between
the look-ahead-bias of using the full data set for imputation and the larger
variance in the imputation from using only the training data. As far as we
know, this paper is the first to introduce a systematic approach to study this
trade-off. We propose a novel Bayesian solution which fuses the use of
globally observed data (i.e. data collected during the full time horizon) and
locally observed data (i.e. data collected up to a fixed time in the past) in
order to optimize the trade-off between look-ahead-bias (which increases when
using globally observed data to impute) and variance (which increases when
using only local data to impute). More specifically, we connect layers of
information revealed in time, through a Bayesian consensus posterior that
fuses an arbitrary number of posteriors to optimally control the variance and
look-ahead-bias trade-off in the imputation.

The use of statistical and machine learning methods in finance has usually
been very careful in avoiding a look-ahead-bias in the downstream task for a
given financial panel without missing values. For example, Gu, Kelly, and Xiu
(2020) and Chen, Pelger, and Zhu (2019) train their machine learning methods
on the training data in the first part of the sample to estimate trading
strategies and evaluate their performance on the out-of-sample test data at
the end of the sample. They follow one common approach in the literature to
use only information from the training data to impute missing information.
Their results could potentially be improved by better imputed values in the
training data if those have a sufficiently small look-ahead-bias. In contrast,
Giglio, Liao, and Xiu (2021) impute a panel of asset returns with matrix
completion methods applied to the full panel of assets returns. Hence, their
subsequent out-of-sample analysis of abnormal asset returns has an implicit
in-sample component as future returns have been used in the imputation. These
two cases of full data imputation (with look-ahead-bias and small variance)
and conservative imputation (without look-ahead- bias and with large variance)
are both special cases of our Bayesian approach. Each of them corresponds to a
different posterior. Our framework allows us to obtain a “consensus” posterior
that is optimal in the sense of this trade-off. Hence, it can be optimal to
use at least some future information if it has a negligible effect on the
look-ahead-bias. This becomes particularly relevant if the amount of missing
data on the training data is large.

As mentioned earlier, we provide the first systematic framework for evaluating
look-ahead-bias and for optimizing the the trade-off between look-ahead-bias
and variance in imputations for future downstream tasks. Any reasonable
formulation of this problem needs to consider three ingredients. First, a
representative downstream task that captures stylized features of broad
interest. Second, a reasonable way to fuse or combine imputation estimators
which are obtained using increasing layers of time information. And, finally,
a way to connect the first and second ingredients to optimally perform the
desired bias-variance trade off. The objective is to find a way to
systematically produce these three ingredients in a tractable and flexible
way. To this end, we adopt a Bayesian imputation framework. This allows us to
precisely define look-ahead-bias as deviation from a specific performance
related to a reasonable downstream task, which we take to be optimizing
expected returns in portfolio selection. The expected return is computed with
the baseline posterior distribution obtained using only locally observed
information. As it is standard in Bayesian imputation (Little and Rubin, 2002;
Rubin, 1987), we impute by sampling (multiple times) from the posterior
distribution given the observed entries by creating multiple versions of
imputed data set which can be used to estimate the optimal posterior return.
Our theoretical results focus on the fundamental case of estimating mean
returns in order to clearly illustrate the novel conceptual ideas, but can be
generalized to more general risk management objectives.

The second ingredient involves a way to fuse inferences arising from different
layers of information. Our Bayesian approach is convenient because we can
optimally combine posterior distributions using a variational approach,
namely, the weighted Barycenter distribution (i.e. the result of an infinite
dimensional optimization problem which finds a consensus distribution which
has the minimal combined discrepancy among a non-parametric family of
distributions based on a given discrepancy criterion). We use the (forward and
backward) Kullback-Leibler divergences and the Wasserstein distance as
criteria to obtain three different ways to generate optimally fused posterior
distributions corresponding to different layers of information (indexed by
time). These optimally fused distributions are indexed by weights
(corresponding to the relative importance given to each distribution to be
combined). It is important to note that since our fusing approach is non-
parametric, this method could be used in combination with other downstream
tasks and is not tied to the particular choice that we used in our first
ingredient.

Finally, the third ingredient in our Bayesian method minimizes the variance of
the optimal return strategy (using the expected return as selection criterion)
assuming a consensus distribution given by the second ingredient - therefore
indexed by the weights attached to each of the posteriors. The optimization is
performed over the such weights under a given look-ahead-bias constraint. The
constrain parameter can be chosen via cross validation or by minimizing the
estimated mean squared error over all consensus posteriors as a function of
the constrain parameter.

In sum, the overall procedure gives the optimal utilization of different
information layers in time, based on a user-defined bias-variance trade-off
relative to the chosen downstream task, and using an optimal consensus
mechanism to fuse the different layers of information.

We illustrate our novel conceptual framework with the data generating
distribution $Z_{t}=\theta+\epsilon_{t}$, where $\theta$ is the mean parameter
constant over time and $\epsilon_{t}$ are the noise terms. For example, in
modeling financial returns, assuming time-stationarity (in particular a
constant mean) over moderate time horizons is a common practice in finance
(Tsay, 2010). The noise $\epsilon_{t}$ can in general follow an $AR(p)$
process. A suitable joint prior on $\theta$ is imposed, after which the joint
posterior conditional on observed entries can in principle be calculated. We
aggregate the different posteriors corresponding to different layers of
information using the variational approach prescribed by a set of weights.
This set of weights is further optimized to find the optimal consensus
posterior according to a bias-variance criterion. This additional step of
optimization over importance weights differentiates our work from the usual
Bayesian fusion setup, as most of the literature only considers a fixed choice
of the weights (Dai, Pollock, and Roberts, 2021; Claici, Yurochkin, Ghosh, and
Solomon, 2020).

To make the two-step optimization procedure tractable, we assume that the
noise $\epsilon_{t}$ are i.i.d Gaussian with a known covariance matrix and use
conjugate prior distributions. This assumption makes the posterior of the
missing entries a conditional Gaussian with mean sampled from the posterior of
$\theta$, and thus we consider only aggregating the different posteriors of
$\theta$. We provide tractable reformulations of the two-step optimization
procedure using the forward KL divergence and Wasserstein distance, when
imposing an additional simultaneous diagonalizability condition on the
covariance matrices of the individual posteriors. While our reformulation
using backward KL divergence does not require simultaneous diagonalizability,
the resulting optimization problem is non-convex and may not be scalable. We
thus advocate using the forward KL divergence and Wasserstein distance over
the backward KL divergence, unless the simultaneous diagonalizability
condition is seriously violated for the type of data at hand.

The contributions of this paper can be summarized as follows:

  1. 1.

We propose a novel conceptual framework to formalize and quantify the look-
ahead-bias in the context of panel data imputation. We combine different
layers of additional future information used in the imputation via optimal
non-parametric consensus mechanisms among Bayesian posterior distributions in
a way that optimally provides a trade-off between the look-ahead-bias and the
variance.

  2. 2.

We illustrate the conceptual framework with a data generating distribution
that has a constant mean vector over time. We propose a two-step optimization
procedure for finding the optimal consensus posterior. The first step of
optimization deals with finding the best possible (non-parametric)
distribution that summarizes the different posterior models given a set of
importance weights. The second step of optimization concerns with finding the
optimal weights, based on a bias-variance criterion that minimizes the
variance while explicitly controlling the degree of the look-ahead-bias.

  3. 3.

To make the two-step procedures tractable, we further assume the noise is
i.i.d multivariate Gaussian with known covariance matrix and use a conjugate
prior distribution. Moreover, we discuss using the (forward and backward) KL-
divergence and Wasserstein distance as the measure of dissimilarity between
posterior distributions. We advise against the backward KL divergence unless
in situations where the simultaneous diagonalizability condition on the
covariance matrices of the individual posteriors is drastically violated.

Finally, we empirically validate the benefits of our methodology with the
downstream task of portfolio optimization. We show that the use of an optimal
consensus posterior results in better performance of the regret of out-of-
sample testing returns measured by what we call “expected conditional mean
squared error”, compared to the use of naive extreme posteriors (corresponding
to focusing either using only the training part of the data for imputation or
the full panel data for imputation). In our companion paper Blanchet,
Hernandez, Nguyen, Pelger, and Zhang (2021), we illustrate the framework to a
special case involving the Wasserstein interpolation of these two extreme
posteriors. We show that the many more layers of additional information to be
fused result in better out-of-sample performance compared to Blanchet,
Hernandez, Nguyen, Pelger, and Zhang (2021). We demonstrate the power of our
method on both simulated and real financial data, and under different
practically relevant missing patterns including missing at random, dependent
block missing, and missing by values.

This paper focuses on the finance setting in which handling missing data is
particularly prevalent, but the distinction between globally and locally
observed data is important in other time-series data as well, for example in
electronic health records. Importantly, it generalizes beyond the time-series
setting to data sets with network structure, where each sample corresponds to
a vertex in the network. The globally observed data is the entire collection
of vertices in the network, while the locally observed data is a subset of the
vertices defined by a neighborhood in this network. One may use the locally
observed data as the training data set for a relevant downstream task while
using the globally observed data for testing. The vertices that are inside of
the neighborhood could potentially have a causal effect (with respect to a
downstream task) on vertices outside, but the causality in the other direction
is not permitted. Such a framework can be used as an approximate model for,
e.g, the prediction of crime rates in different parts of a state, of the
radiation level of surrounding cities after a damage of nuclear power plant,
or of the earthquake magnitude in a geographic region. If missing data is
imputed globally using all information, an out-of-sample performance of the
downstream task will inherently suffer from overfitting as the imputation
extracts undesirable information. The imputation using only locally observed
data again fails to utilize all valuable information to improve the variance.

This paper is organized as follows. In Section 2 we discuss the related
literature on imputation methods. In Section 3 we extend the bias-variance
multiple imputation framework to more than two individual posteriors along
with two basic modules of the framework: the posterior generator and the
multiple imputation sampler. Section 4 discusses the consensus mechanism
module which complements our framework. Subsequently, in Section 5 and 7, we
use the forward and the backward KL divergence as a measure of dissimilarity
between distributions and derive a reformulation of the consensus mechanism,
while in Section 6 we use the Wasserstein distance. Section 8 reports the
simulation and empirical results. The Appendix collects the proofs and
additional empirical results.

##  2 Related Literature

Strategies for missing data imputation can be classified into two categories:
single and multiple imputation. Single imputation refers to the generation and
usage of a specific number (i.e., a best guess) in place of each missing
value. Standard non-parametric methods for single imputation include
$k$-nearest neighbor, mean/median imputation, smoothing, interpolation and
splines (Kreindler and Lumsden, 2006). Alternatively, single imputation can
also be achieved using a parametric approach such as the popular EM algorithm
of Dempster, Laird, and Rubin (1977), Bayesian networks (Lan, Xu, Ma, and Li,
2020) or $t$-copulas and transformed $t$-mixture models (Friedman, Huang,
Zhang, and Cao, 2012). Single imputation also includes the use of iterative
methods in multivariate statistics literature to deal with missing values
encountered in classical tasks such as correspondence analysis (de Leeuw and
van der Heijden (1988)), multiple correspondence analysis, and multivariate
analysis of mixed data sets (Audigier, Husson, and Josse, 2016). Recently,
there is a surge in the application of recurrent neural networks for
imputation, as well as generative adversarial networks.111Examples of
recurrent networks for imputation include Che, Purushotham, Cho, Sontag, and
Liu (2018); Choi, Bahadori, Schuetz, Stewart, and Sun (2016); Lipton, Kale,
and Wetzel (2016); Cao, Wang, Li, Zhou, Li, and Li (2018), while Yoon, Jordon,
and van der Schaar (2018) is an example for the GAN approach.

A closely related field of single imputation of panel data is matrix
completion and matrix factorization.222This includes the work of Cai, Candès,
and Shen (2010); Candès and Recht (2009); Candès and Tao (2010); Mazumder,
Hastie, and Tibshirani (2010); Keshavan, Montanari, and Oh (2010). Chen, Fan,
Ma, and Yan (2019) propose a de-biased estimator and an optimal construction
of confidence intervals/regions for the missing entries and the low-rank
factors. Xiong and Pelger (2019) impute missing data with a latent factor
model by applying principal component analysis to an adjusted covariance
matrix estimated from partially observed panel data. They provide an
asymptotic inferential theory and allow for general missing patterns, which
can also depend on the underlying factor structure. This is important as
Bryzgalova, Lettau, Lerner, and Pelger (2022) show in their comprehensive
empirical analysis, that missing patterns of firm specific and asset pricing
information are impacted by the cross- and time-series dependency structure in
the data. Bai and Ng (2020) suggest an imputation procedure that uses factors
estimated from a tall block along with the re-rotated loadings estimated from
a wide block.

Multiple imputation generalizes the single imputation procedures in that the
missing entries are filled-in with multiple guesses instead of one single
guess or estimate, accounting for the uncertainty involved in imputation. The
use of multiple imputation results in $m$ ($m>2$) completed data sets, each is
then processed and analyzed with whatever statistical method can be used as if
there is no missing data. Rubin’s rule (Rubin, 1987) is often used to combine
the quantity of interest computed from each completed data set, which averages
over the point estimates and uses a slightly more involved expression for the
standard errors.

The central idea underpinning the multiple imputation is the Bayesian
framework. This framework samples the missing values from a posterior
predictive distribution of the missing entries given the observed data.
Multiple imputation is first developed for non-responses in survey sampling
(Rubin, 1987). Since then, it has been extended to time-series data (Honaker
and King, 2010; Halpin, 2016), where a key new element is to preserve
longitudinal consistency in imputation. More recently, Honaker and King (2010)
uses smooth basis functions to increase smoothness of the imputation values
while Halpin (2016) uses gap-filling sequence of imputation for categorical
time-series data and produces smooth patterns of transitions.333The
flexibility and predictive performance of multiple imputation method have been
successfully demonstrated in a variety of data sets, such as Industry and
Occupation Codes (Clogg, Rubin, Schenker, Schultz, and Weidman, 1991), GDP in
African nation (Honaker and King, 2010), and concentrations of pollutants in
the arctic (Hopke, Liu, and Rubin, 2001).

This paper differs from standard multiple imputation literature in that we are
interested in achieving an optimal trade-off between look-ahead-bias and
variance by unifying the posteriors corresponding to different layers of
information into a single “consensus” posterior, while most of the literature
in multiple imputation would advocate using the full panel data. The
unification of different posteriors into a single coherent inference, or
(Bayesian) “fusion”, has recently gained considerable interest in the research
community (Dai, Pollock, and Roberts, 2021; Claici, Yurochkin, Ghosh, and
Solomon, 2020), but often the different posteriors arise from the context of
using paralleled computing machines (Scott, Blocker, Bonassi, Chipman, George,
and McCulloch, 2016), multiple sensors in multitarget tracking (Li, Hlawatsch,
and Djurić, 2019; Da, Li, Zhu, Fan, and Fu, 2020), or in model selection
(Buchholz, Ahfock, and Richardson, 2019). In these settings, typically all of
the distributions are given equal weight. Our approach contributes to the
fusion literature by considering a downstream criterion that can be used to
optimize the importance (or weight) of each distribution.

##  3 The Bias-Variance Multiple Imputation Framework

Notations. For any integers $K\leq L$, we use $[K]$ to denote the set
$\\{1,\ldots,K\\}$ and $[K,L]$ to denote the set $\\{K,\ldots,L\\}$. We use
$\mathcal{M}$ to denote the space of all probability measures supported on
$\mathbb{R}^{n}$. The set of all $n$-dimensional Gaussian measures is denoted
by $\mathcal{N}_{n}$, and we use $\mathbb{S}_{++}^{n}$ to represent the set of
positive-definite matrices of dimension $n\times n$.

We assume a universe of $n$ assets over $T$ periods. The joint return of these
$n$ assets at any specific time $t$ is represented by a random vector
$Z_{t}\in\mathbb{R}^{n}$. The random vector $M_{t}\in\\{0,1\\}^{n}$ indicates
the pattern of the missing values at time $t$, that is,

$\forall i\in[n]:\qquad(M_{t})_{i}=\begin{cases}1&\text{if the $i$-th
component of $Z_{t}$ is missing},\\\ 0&\text{if the $i$-th component of
$Z_{t}$ is observed}.\end{cases}$

Coupled with the indicator random variables $M_{t}$, we define the following
projection operator $\mathcal{P}_{M_{t}}$ that projects444 The projection
operator $\mathcal{P}_{M_{t}}$ can be defined using a matrix-vector
multiplication as follows. Upon observing $M_{t}$, let
$\mathcal{I}=\\{i\in[n]:(M_{t})_{i}=1\\}$ be the set of indices of the missing
entries. We assume without any loss of generality that $\mathcal{I}$ is
ordered and has $d$ elements, and that the $j$-th element of $\mathcal{I}$ can
be assessed using $\mathcal{I}_{j}$. Define the matrix
$P_{M_{t}}\in\\{0,1\\}^{d\times n}$ with its $(j,\mathcal{I}_{j})$ element
being 1 for $j\in[d]$, and the remaining elements are zeros, then
$\mathcal{P}_{M_{t}}(Z_{t})=P_{M_{t}}Z_{t}$. The operator
$\mathcal{P}_{M_{t}}^{\perp}$ can be defined with a matrix $P_{M_{t}}^{\perp}$
in a similar manner using the set
$\mathcal{I}^{\perp}=\\{i\in[n]:(M_{t})_{i}=0\\}$. the latent variables
$Z_{t}$ to its unobserved components $Y_{t}$

$\mathcal{P}_{M_{t}}(Z_{t})=Y_{t}.$

Consequentially, the orthogonal projection $\mathcal{P}_{M_{t}}^{\perp}$
projects $Z_{t}$ to the observed components $X_{t}$

$\mathcal{P}_{M_{t}}^{\perp}(Z_{t})=X_{t}.$

At the fundamental level, we assume the following generative model

$\forall t:\quad\left\\{\begin{array}[]{l}Z_{t}=\theta+\varepsilon_{t}\\\
X_{t}=\mathcal{P}_{M_{t}}^{\perp}(Z_{t}),~{}Y_{t}=\mathcal{P}_{M_{t}}(Z_{t}),\end{array}\right.\qquad\text{where}\qquad\varepsilon_{t}\in\mathbb{R}^{n},~{}M_{t}\in\\{0,1\\}^{n}.$
(1)

When the generative model (1) is executed over the $T$ periods, the observed
variables $(X_{t},M_{t})$ are accumulated to produce a data set $\mathcal{D}$
of $n$ rows and $T$ columns. This data set is separated into a training set
consisting of $T^{train}$ periods, and the remaining $T^{test}=T-T^{train}$
periods are used as the testing set. We will be primarily concerned with the
imputation of the training set.

Figure 1: Illustration of panel with missing data This figure illustrates
panel data $\mathcal{D}$ with missing entries in black. The observed
components of $Z_{t}$ are denoted by $X_{t}$, and the unobserved components of
$Z_{t}$ are denoted by $Y_{t}$.

In this paper, we pursue the Bayesian approach: we assume that the vector
$\theta\in\mathbb{R}^{n}$ is unknown, and is treated as a Bayesian parameter
with a prior distribution $\pi_{0}$; we also treat $(Y_{t})_{t\in[T]}$ as
Bayesian parameters with appropriate priors. Based on the data set with
missing entries $\mathcal{D}$, our strategy relies on computing the
distribution of $(Y_{t})_{t\in[T]}$ conditional on the data set $\mathcal{D}$,
and then generate multiple imputations of $(Y_{t})_{t\in[T]}$ by sampling from
this distribution. To compute the posterior of
$(Y_{t})_{t\in[T]}|\mathcal{D}$, we first need to calculate a posterior
distribution of the unknown mean vector $\theta$ conditional on $\mathcal{D}$,
and then the distribution of $(Y_{t})_{t\in[T]}|\mathcal{D}$ is obtained by
marginalizing out $\theta$ from the distribution of
$(Y_{t})_{t\in[T]}|(\theta,\mathcal{D})$.

We endeavor in this paper to explore a novel approach to generate multiply-
imputed version of the data set. We make the following assumptions.

######  Assumption 3.1 (Observability).

We observe at least one value for each row in the training section of
$\mathcal{D}$.

######  Assumption 3.2 (Ignorability).

The missing pattern satisfies the ignorability assumption (Little and Rubin,
2002), namely, the probability of $(M_{t})_{i}=1~{}\forall i\in[n]~{}\forall
t\in[T]$ does not depend on $(Y_{t})_{t\in[T]}$ or $\theta$, conditional on
$(X_{t})_{t\in[T]}$.

######  Assumption 3.3 (Bayesian conjugacy).

The noise $\varepsilon_{t}~{}\forall t\in[T]$ in the latent generative model
(1) is independently and identically distributed as
$\mathcal{N}_{n}(0,\Omega)$, where $\Omega\in\mathbb{S}_{++}^{n}$ is a known
covariance matrix. The prior $\pi_{0}$ is either a non-informative flat prior
or Gaussian with mean vector $\mu_{0}\in\mathbb{R}^{n}$ and covariance matrix
$\Sigma_{0}\in\mathbb{S}_{++}^{n}$. The priors on $Y_{t}$ conditional on
$\theta$ are independent across $t$, and the conditional distribution is
Gaussian $\mathcal{N}_{\mathrm{dim}(Y_{t})}(\theta_{Y_{t}},\Omega_{Y_{t}})$,
where $\theta_{Y_{t}}=\mathcal{P}_{M_{t}}(\theta)$, and $\Omega_{Y_{t}}$ is
obtained by applying the projection operator $\mathcal{P}_{M_{t}}$ on
$\Omega$.555Using the construction of $P_{M_{t}}$ in Footnote 4, we have
$\Omega_{Y_{t}}=P_{M_{t}}\Omega P_{M_{t}}^{\top}$. Likewise for
$\Omega_{X_{t}}$. The latent generative model (1) implies that the likelihood
of $X_{t}$ conditional on $(\theta,Y_{t})$ is Gaussian
$N_{\mathrm{dim}(X_{t})}\big{(}\theta_{X_{t}}+\Omega_{X_{t}Y_{t}}\Omega_{Y_{t}}^{-1}(Y_{t}-\theta_{Y_{t}}),\Omega_{X_{t}}-\Omega_{X_{t}Y_{t}}\Omega_{Y_{t}}^{-1}\Omega_{Y_{t}X_{t}}\big{)}$,
where $\theta_{X_{t}}=\mathcal{P}_{M_{t}}^{\perp}(\theta)$, the matrix
$\Omega_{Y_{t}X_{t}}\in\mathbb{R}^{\mathrm{dim}(Y_{t})\times\mathrm{dim}(X_{t})}$
is the covariance matrix between $Y_{t}$ and $X_{t}$ induced by $\Omega$,
similarly for $\Omega_{X_{t}Y_{t}}$.

Following our discussion on the look-ahead bias, we can observe that the
distribution of $\theta|\mathcal{D}$ may carry look-ahead bias originating
from the testing portion of $\mathcal{D}$, and this look-ahead bias will be
internalize into the look-ahead bias of
$(Y_{t})_{t\in[T^{train}]}|\mathcal{D}$. To mitigate this negative effect, it
is crucial to re-calibrate the distribution of $\theta|\mathcal{D}$ to reduce
the look-ahead bias on the posterior distribution of the mean parameter
$\theta$, with the hope that this mitigation will be spilt-over to a similar
mitigation of the look-ahead bias on the imputation of
$(Y_{t})_{t\in[T^{train}]}$. Our BVMI framework pictured in Figure 2 is
designed to mitigate the look-ahead bias by implementing there modules:

  * (G)

The posterior Generator

  * (C)

The Consensus mechanism

  * (S)

The multiple imputation Sampler

Figure 2: Illustration of workflow This figures illustrates the workflow of
our framework. The Bias-Variance Multiple Imputation (BVMI) framework
(contained in dashed box) receives a data set with missing values as an input,
and generates multiple data sets as outputs. The posterior generator (G) can
exploit time dependency, while the consensus mechanism (C) is bias-variance
targeted.

For the rest of this section, we will elaborate on the generator (G) in
Subsection 3.1, and we detail the sampler (S) in Subsection 3.2. The
construction of the consensus mechanism (C) is more technical, and thus the
full details on (C) will be provided in subsequent sections.

###  3.1 Posterior Generator

We now discuss in more detail how to generate multiple posterior distributions
by taking into account the time-dimension of the data set. To this end, we fix
the desired number of posterior distributions to $K$, ($2\leq K\leq
T^{test}+1)$. Moreover, we assign the integer values $T_{1},\ldots,T_{K}$
satisfying

$T^{train}=T_{1}<T_{2}<\ldots<T_{K}=T$

to indicate the amount of data that is synthesized in each posterior. The data
set truncated to time $T_{k}$ is denoted as
$\mathcal{D}_{k}=\\{(X_{t},M_{t}),t\in[T_{k}]\\}$, note that $\mathcal{D}_{1}$
coincides with the training data set while $\mathcal{D}_{K}$ coincides with
the whole data set $\mathcal{D}$. By construction, the sets $\mathcal{D}_{k}$
inherits a hierarchy in terms of time: for any $k<k^{\prime}$, the data set
$\mathcal{D}_{k}$ is part of $\mathcal{D}_{k^{\prime}}$, and thus
$\mathcal{D}_{k^{\prime}}$ contains at least as much of information as
$\mathcal{D}_{k}$. The $k$-th posterior for $\theta$ is formulated conditional
on $\mathcal{D}_{k}$. Under standard regularity conditions, the posterior
$\pi_{k}$ can be computed from the prior distribution $\pi_{0}$ using the
Bayes’ theorem (Schervish, 1995, Theorem 1.31) as

$\frac{\mathrm{d}\pi_{k}}{\mathrm{d}\pi_{0}}(\theta|\mathcal{D}_{k})=\frac{f(\mathcal{D}_{k}|\theta)}{\int_{\mathbb{R}^{n}}f(\mathcal{D}_{k}|\vartheta)\pi_{0}(\mathrm{d}\vartheta)},$
(2)

where $\mathrm{d}\pi_{k}/\mathrm{d}\pi_{0}$ represents the Radon-Nikodym
derivative of $\pi_{k}$ with respect to $\pi_{0}$, and
$f(\mathcal{D}_{k}|\theta)$ is the conditional density representing the
likelihood of observing $\mathcal{D}_{k}$ given the parameter $\theta$. Under
the generative model (1), and Assumption 3.2, we have

$f(\mathcal{D}_{k}|\theta)=f_{X}\left((X_{t})_{t\in[T_{k}]}|\theta\right)\times
f_{M}\left((M_{t})_{t\in[T_{k}]}|(X_{t})_{t\in[T_{k}]},\theta\right)=f_{X}\left((X_{t})_{t\in[T_{k}]}|\theta\right)\times\tilde{f}_{M}\left((M_{t})_{t\in[T_{k}]}|(X_{t})_{t\in[T_{k}]}\right)$

for some appropriate (conditional) likelihood functions $f_{X},f_{M}$ and
$\tilde{f}_{M}$. In this case, the terms involving $(M_{t})_{t\in[T_{k}]}$ in
both the numerator and denominator of equation (2) cancel each other out,
hence (2) can be simplified to

$\frac{\mathrm{d}\pi_{k}}{\mathrm{d}\pi_{0}}(\theta|\mathcal{D}_{k})=\frac{f_{X}\left((X_{t})_{t\in[T_{k}]}|\theta\right)}{\int_{\mathbb{R}^{n}}f_{X}\left((X_{t})_{t\in[T_{k}]}|\vartheta\right)\pi_{0}(\mathrm{d}\vartheta)},$

Under Assumption 3.3, the posterior distribution $\pi_{k}$ can be computed
more conveniently. Assumption 3.3 implies that conditioning on $\theta$, the
latent variables $Z_{t}$ and $Z_{t^{\prime}}$ are i.i.d. following a Gaussian
distribution with mean $\theta$ and covariance matrix $\Omega$ for any $t$,
$t^{\prime}$. Because the operator $\mathcal{P}_{M_{t}}^{\perp}$ is an affine
transformation (see Footnote 4), Theorem 2.16 in Fang, Kotz, and Ng (1990)
asserts that $X_{t}|\theta$ is Gaussian with mean $\theta_{X_{t}}$ and
covariance matrix $\Omega_{X_{t}}$. In this case, the posterior $\pi_{k}$ can
be computed as

$\frac{\mathrm{d}\pi_{k}}{\mathrm{d}\pi_{0}}(\theta|\mathcal{D}_{k})=\frac{\prod_{t\in[T_{k}]}\mathcal{N}_{\mathrm{dim}(X_{t})}(X_{t}|\theta_{X_{t}},\Omega_{X_{t}})}{\int_{\mathbb{R}^{n}}\prod_{t\in[T_{k}]}\mathcal{N}_{\mathrm{dim}(X_{t})}(X_{t}|\vartheta_{X_{t}},\Omega_{X_{t}})\pi_{0}(\mathrm{d}\vartheta)},$

where $\mathcal{N}_{\mathrm{dim}(X_{t})}(\cdot|\theta_{X_{t}},\Omega_{X_{t}})$
is the probability density function of the Gaussian distribution
$\mathcal{N}_{\mathrm{dim}(X_{t})}(\theta_{X_{t}},\Omega_{X_{t}})$. Because
$\pi_{0}$ is either flat or Gaussian by Assumption 3.3, the Bayes’ theorem
asserts that the posterior $\pi_{k}$ also admits a density (with respect to
the Lebesgue measure on $\mathbb{R}^{n}$). This density (we first take the
case that $\pi_{0}$ is Gaussian for illustration) evaluated at any point
$\theta\in\mathbb{R}^{n}$ is proportional to

$\displaystyle\mathcal{N}_{n}(\theta|\mu_{0},\Sigma_{0})\prod_{t\in[T_{k}]}\mathcal{N}_{\mathrm{dim}(X_{t})}(X_{t}|\theta_{X_{t}},\Omega_{X_{t}})$
$\displaystyle\propto$
$\displaystyle\exp\left(-\frac{1}{2}(\theta-\mu_{0})^{\top}\Sigma_{0}^{-1}(\theta-\mu_{0})\right)\prod_{t\in[T_{k}]}\exp\left(-\frac{1}{2}(\theta_{X_{t}}-X_{t})^{\top}\Omega_{X_{t}}^{-1}(\theta_{X_{t}}-X_{t})\right).$

By expanding the exponential terms and completing the squares, the posterior
distribution $\theta|\mathcal{D}_{k}$ is then governed by
$\pi_{k}\sim\mathcal{N}_{n}(\mu_{k},\Sigma_{k})$ with

$\Sigma_{k}=\left(\Sigma_{0}^{-1}+\sum_{t\in[T_{k}]}(\mathcal{P}_{M_{t}}^{\perp})^{-1}(\Omega_{X_{t}}^{-1})\right)^{-1},\quad\mu_{k}=\Sigma_{k}\left(\Sigma_{0}^{-1}\mu_{0}+\sum_{t\in[T_{k}]}(\mathcal{P}_{M_{t}}^{\perp})^{-1}(\Omega_{X_{t}}^{-1})\
(\mathcal{P}_{M_{t}}^{\perp})^{-1}(X_{t})\right),$ (3)

where $(\mathcal{P}_{M_{t}}^{\perp})^{-1}$ is the inverse map of
$\mathcal{P}_{M_{t}}^{\perp}$ obtained by filling entries with zeros666This
inverse map can be defined as
$(\mathcal{P}_{M_{t}}^{\perp})^{-1}(\Omega_{X_{t}}^{-1})=\arg\min\\{\|\tilde{\Omega}\|_{0}:\tilde{\Omega}\in\mathbb{R}^{n\times
n},~{}\mathcal{P}_{M_{t}}^{\perp}(\tilde{\Omega})=\Omega_{X_{t}}^{-1}\\}$,
where $\|\cdot\|_{0}$ is a sparsity inducing norm. Similarly for
$(\mathcal{P}_{M_{t}}^{\perp})^{-1}(X_{t})$.. An elementary linear algebra
argument asserts that $\Omega_{X_{t}}^{-1}$ is positive definite, and hence
$(\mathcal{P}_{M_{t}}^{\perp})^{-1}(\Omega_{X_{t}}^{-1})$ is positive
semidefinite. The matrix inversion that defines $\Sigma_{k}$ is thus well-
defined thanks to the positive definiteness of $\Sigma_{0}$ in Assumption 3.3.
Alternatively, in Assumption 3.3 if we impose a flat uninformative prior on
$\theta$, then the posterior distribution of $\theta|\mathcal{D}_{k}$ is still
Gaussian with $\pi_{k}\sim\mathcal{N}_{n}(\mu_{k},\Sigma_{k})$ where

$\Sigma_{k}=\left(\sum_{t\in[T_{k}]}(\mathcal{P}_{M_{t}}^{\perp})^{-1}(\Omega_{X_{t}}^{-1})\right)^{-1},\quad\mu_{k}=\Sigma_{k}\left(\sum_{t\in[T_{k}]}(\mathcal{P}_{M_{t}}^{\perp})^{-1}(\Omega_{X_{t}}^{-1})\
(\mathcal{P}_{M_{t}}^{\perp})^{-1}(X_{t})\right).$ (4)

The following lemma ensures us that the covariance $\Sigma_{k}$ is positive
definite.

######  Lemma 3.4.

Under Assumption 3.1, the posterior covariance (4) is positive definite.

The product of the posterior generator module $(G)$ is $K$ elementary
posteriors $\\{\pi_{k}\\}_{k\in[K]}$, where $\pi_{k}$ is the posterior of
$\theta|\mathcal{D}_{k}$. For any $k<k^{\prime}$, the definition of the set
$\mathcal{D}_{k}$ and $\mathcal{D}_{k^{\prime}}$ implies that the posterior
$\pi_{k^{\prime}}$ carries at least as much of look-ahead bias as the
posterior $\pi_{k}$. The collection of posteriors $\\{\pi_{k}\\}$ is
transmitted to the consensus mechanism $(C)$ to form an aggregated posterior
$\pi^{\star}$ that strikes a balance between the variance and the look-ahead-
bias. The aggregated posterior $\pi^{\star}$ is injected to the sampler module
(S), which we now study.

###  3.2 Multiple Imputation Sampler

A natural strategy to impute the missing values depends on recovering the
(joint) posterior distribution of $(Y_{t})_{t\in[T]}$ conditional on the
observed data $\mathcal{D}$. Fortunately, if the noise $\epsilon_{t}$ are
mutually independent, and the priors on $Y_{t}$ are mutually independent given
$\theta$, it suffices to consider the posterior distribution of $Y_{t}$ for
any given $t$. Using the law of conditional probability and leveraging the
availability of the aggregated posterior $\pi^{\star}$, we have for any
measurable set $\mathcal{Y}\subseteq\mathbb{R}^{n}$

$\displaystyle\mathbb{P}(Y_{t}\in\mathcal{Y}|\mathcal{D})$
$\displaystyle=\int_{\mathbb{R}^{n}}\mathbb{P}(Y_{t}\in\mathcal{Y}|\theta,\mathcal{D})\pi^{\star}(\mathrm{d}\theta)$
$\displaystyle=\int_{\mathbb{R}^{n}}\mathbb{P}(Y_{t}\in\mathcal{Y}|\theta,X_{t},M_{t})\pi^{\star}(\mathrm{d}\theta),$
(5a) where the last equality holds thanks to the independence between $Y_{t}$
and any $(X_{t^{\prime}},M_{t^{\prime}})$ for $t^{\prime}\neq t$ conditional
on $(\theta,X_{t},M_{t})$. Evaluating the posterior distribution for
$Y_{t}|\mathcal{D}$ thus entails a multi-dimensional integration, which is in
general computationally intractable.

Under Assumption 3.3, $Z_{t}$ conditional on $\theta$ is Gaussian
$\mathcal{N}_{n}(\theta,\Omega)$. Because $Z_{t}$ is decomposed into $X_{t}$
and $Y_{t}$, the distribution of $Y_{t}$ conditional on $\theta$ and
$(X_{t},M_{t})$ is thus Gaussian (Cambanis, Huang, and Simons, 1981, Corollary
5). More specifically,

$Y_{t}|\theta,X_{t},M_{t}\sim\mathcal{N}_{\mathrm{dim}(Y_{t})}\left(\theta_{Y_{t}}+\Omega_{Y_{t}X_{t}}\Omega_{X_{t}}^{-1}(X_{t}-\theta_{X_{t}}),\Omega_{Y_{t}}-\Omega_{Y_{t}X_{t}}\Omega_{X_{t}}^{-1}\Omega_{X_{t}Y_{t}}\right),$
(5b)

where $\theta_{X_{t}}=\mathcal{P}_{M_{t}}(\theta)$,
$\theta_{Y_{t}}=\mathcal{P}_{M_{t}}^{\perp}(\theta)$, while $\Omega_{X_{t}}$
and $\Omega_{Y_{t}}$ are obtained by applying the projection operators
$\mathcal{P}_{M_{t}}$ and $\mathcal{P}_{M_{t}}^{\perp}$ on $\Omega$. The
matrix
$\Omega_{Y_{t}X_{t}}\in\mathbb{R}^{\mathrm{dim}(Y_{t})\times\mathrm{dim}(X_{t})}$
is the covariance matrix between $Y_{t}$ and $X_{t}$ induced by $\Omega$.
Notice that by the definitions of the projection operators, $\theta_{X_{t}}$
and $\theta_{Y_{t}}$ are affine transformations of $\theta$. As a consequence,
the mean of $Y|\theta,X_{t},M_{t}$ is an affine function of $\theta$, while
its covariance matrix is independent of $\theta$. The posterior formula (5a)
coupled with the structural form (5b) indicate that the distribution of
$Y_{t}|\mathcal{D}$ coincide with the distribution of the random vector
$\xi_{t}\in\mathbb{R}^{\mathrm{dim}(Y_{t})}$ dictated by

$\xi_{t}=A_{t}\theta+b_{t}+\eta_{t},\qquad\theta\sim\pi^{\star},\qquad\eta_{t}\sim\mathcal{N}_{\mathrm{dim}(Y_{t})}(0,S_{t}),\qquad\theta\perp\\!\\!\\!\perp\eta_{t},$
(6)

for some appropriate parameters $A_{t}$, $b_{t}$ and $S_{t}$ that are designed
to match the mean vector and the covariance matrix of the Gaussian
distribution in (5b). More specifically, we have
$A_{t}\theta=\theta_{Y_{t}}-\Omega_{Y_{t}X_{t}}\Omega_{X_{t}}^{-1}\theta_{X_{t}},b_{t}=\Omega_{Y_{t}X_{t}}\Omega_{X_{t}}^{-1}X_{t},S_{t}=\Omega_{Y_{t}}-\Omega_{Y_{t}X_{t}}\Omega_{X_{t}}^{-1}\Omega_{X_{t}Y_{t}}$.
Algorithm 1 details the procedure for missing value imputation using the model
(6).

The variability in the sampling of $\xi_{t}$ using (6) comes from two sources:
that of sampling $\theta$ from $\pi^{\star}$ and that of sampling $\eta_{t}$
from $\mathcal{N}_{dim(Y_{t})}(0,S_{t})$. As we have seen from Section 3.1
(e.g., equations (3) and (4)), the posterior covariance of $\theta$ is
inverse-proportional to the time-dimension of the data set, and therefore the
variability due to $\pi^{\star}$ is likely to be overwhelmed by the
variability due to $\eta_{t}$ unless $T$ is small. To fully exhibit the
potential of the aggregated posterior $\pi^{\star}$ for bias and variance
trade-off, especially when validating our proposed method in numerical
experiments, we choose to eliminate the idiosyncratic noise $\eta_{t}$,
namely, the missing values are imputed by its conditional expectation given
$\theta$, which is depicted in Algorithm 2.

Algorithm 1 Fully Bayesian Imputation

0: aggregated posterior $\pi^{\star}$, data set
$\mathcal{D}=\\{(X_{t},M_{t}):t\in[T]\\}$, covariance matrix $\Omega$

Sample $\theta$ from $\pi^{\star}$

for $t=1,\ldots,T$ do

Compute $A_{t}$, $b_{t}$, $S_{t}$ for model (6) dependent on
$(X_{t},M_{t},\Omega)$

Sample $\eta_{t}$ from $\mathcal{N}_{\mathrm{dim}(Y_{t})}(0,S_{t})$

Impute $Y_{t}$ by $A_{t}\theta+b_{t}+\eta_{t}$

end for

An imputed data set

Algorithm 2 Conditional Expectation Bayesian Imputation

0: aggregated posterior $\pi^{\star}$, data set
$\mathcal{D}=\\{(X_{t},M_{t}):t\in[T]\\}$, covariance matrix $\Omega$

Sample $\theta$ from $\pi^{\star}$

for $t=1,\ldots,T$ do

Compute $A_{t}$, $b_{t}$ for model (6) dependent on $(X_{t},M_{t},\Omega)$

Impute $Y_{t}$ by $A_{t}\theta+b_{t}$

end for

An imputed data set

##  4 Bias-Variance Targeted Consensus Mechanism

We devote this section to provide the general technical details on the
construction of the consensus mechanism module (C) in our BVMI framework.
Conceptually, the consensus mechanism module receives $K$ posteriors
$\\{\pi_{k}\\}_{k\in[K]}$ from the posterior generator module (G), then
synthesizes a unique posterior $\pi^{\star}$ and transmits $\pi^{\star}$ to
the sampler module (S). We first provide a formal definition of a consensus
mechanism.

######  Definition 4.1 (Consensus mechanism).

Fix an integer $K\geq 2$. A consensus mechanism
$\mathcal{C}_{\lambda}:\mathcal{M}^{K}\to\mathcal{M}$ parametrized by
$\lambda$ outputs a unifying posterior $\pi^{\star}$ for any collection of $K$
elementary posteriors $\\{\pi_{k}\\}$.

We assume that the parameter $\lambda$ is finite dimensional, and the set of
all feasible parameters is denoted by $\Delta_{K}$, which is a nonempty finite
dimensional set. The space of all possible consensus mechanism induced by
$\Delta_{K}$ is represented by
$\mathbb{C}\triangleq\\{\mathcal{C}_{\lambda}:\lambda\in\Delta_{K}\\}$. For a
given set of posteriors $\\{\pi_{k}\\}$, we need to identify the optimal
consensus mechanism $\mathcal{C}_{\lambda^{\star}}\in\mathbb{C}$ in the sense
that the aggregated posterior
$\pi^{\star}=\mathcal{C}_{\lambda^{\star}}(\\{\pi_{k}\\})$ achieves certain
optimality criteria. In this paper, the optimality criteria of $\pi^{\star}$
is determined based on the notion of bias-variance trade-off, which is
ubiquitously used in (Bayesian) statistics. The first step towards this goal
entails quantifying the bias of an aggregated posterior $\pi^{\star}$.

######  Definition 4.2 ($(\ell,\mu,\delta)$-bias tolerable posterior).

Fix a penalization function $\ell:\mathbb{R}^{n}\to\mathbb{R}$, an anchored
mean $\mu\in\mathbb{R}^{n}$ and a tolerance $\delta\in\mathbb{R}_{+}$. A
posterior $\pi^{\star}$ is said to be $(\ell,\mu,\delta)$-bias tolerable if
$\ell(\mathds{E}_{\pi^{\star}}[\theta]-\mu)\leq\delta$.

If $\ell$ is the Euclidean norm and the anchored mean is set to
$\theta^{true}$, then we recover the standard measure of bias in the
literature with the (subjective) posterior distribution in place of the
(objective) frequentist population distribution, and a posterior $\pi^{\star}$
is unbiased if it is $(\|\cdot\|_{2},\theta^{true},0)$-bias tolerable.

Definition 4.2 provides a wide spectrum of options to measure the bias of the
aggregated posterior. Most importantly, we can flexibly choose a suitable
value for the anchored mean $\mu$ for our purpose of mitigating look-ahead
bias. The tolerance parameter $\delta$ indicates how much of look-ahead bias
is accepted in the aggregated posterior $\pi^{\star}$. If $\delta=0$, then we
are aiming for an aggregated posterior distribution that has the exact mean
vector $\mu$. The generality of Definition 4.2 also allows us to choose the
penalization function. Throughout this paper, we will focus on the
penalization function of the form

$\ell_{\mathcal{Z}}(\mu)=\sup\limits_{z\in\mathcal{Z}}~{}z^{\top}\mu$ (7)

parametrized by some nonempty, compact set
$\mathcal{Z}\subseteq\mathbb{R}^{n}$. The functional form (7) enables us to
model diverse bias penalization effects as demonstrated in the following
example.

######  Example 4.3 (Forms of $\ell_{\mathcal{Z}}$).

  1. (i)

If $\mathcal{Z}=\\{z\in\mathbb{R}^{n}:\|z\|_{2}\leq 1\\}$ is the Euclidean
ball, then $\ell_{\mathcal{Z}}(\mu)=\|\mu\|_{2}$. This 2-norm penalization
coincides with the canonical bias measurement in the literature.

  2. (ii)

If $\mathcal{Z}$ is a convex polyhedron
$\mathcal{Z}=\\{z\in\mathbb{R}^{n}:Az\leq b\\}$ for some matrix $A$ and vector
$b$ of appropriate dimensions, then $\ell_{\mathcal{Z}}$ coincides with the
support function of $\mathcal{Z}$.

  3. (iii)

If $\mathcal{Z}=\\{\mathbbm{1}_{n}\\}$, then $\ell_{\mathcal{Z}}$ is an upward
bias penalization under the $1/n$-uniform portfolio. The $1/n$-uniform
portfolio has been shown to be a robust portfolio (DeMiguel, Garlappi, and
Uppal, 2007).

Next, we discuss the variance properties of the aggregated posterior. Let
$\mathrm{Cov}_{\mathcal{C}_{\lambda}(\\{\pi_{k}\\})}[\theta]\in\mathbb{R}^{n\times
n}$ be the covariance matrix of $\theta$ under the posterior distribution
$\mathcal{C}_{\lambda}(\\{\pi_{k}\\})$. We consider the following minimal
variance criterion which is computed based on taking the sum of the diagonal
entries.

######  Definition 4.4 (Minimal variance posterior).

For a fixed input $\\{\pi_{k}\\}$, the aggregated posterior $\pi^{\star}$
obtained by $\pi^{\star}=\mathcal{C}_{\lambda^{\star}}(\\{\pi_{k}\\})$ for an
$\mathcal{C}_{\lambda^{\star}}\in\mathbb{C}$ has minimal variance if

$\mathrm{Tr}(\mathrm{Cov}_{\pi^{\star}}[\theta])\leq\mathrm{Tr}(\mathrm{Cov}_{\mathcal{C}_{\lambda}(\\{\pi_{k}\\})}[\theta])\quad\forall\mathcal{C}_{\lambda}\in\mathbb{C}.$

We are now ready to define our notion of optimal consensus mechanism.

######  Definition 4.5 (Optimal consensus mechanism).

Fix an input $\\{\pi_{k}\\}$. The consensus mechanism
$\mathcal{C}_{\lambda^{\star}}\in\mathbb{C}$ is optimal for $\\{\pi_{k}\\}$ if
the aggregated posterior
$\pi^{\star}=\mathcal{C}_{\lambda^{\star}}(\\{\pi_{k}\\})$ is
$(\ell,\mu,\delta)$-bias tolerable and has minimal variance.

Intuitively, the optimal consensus mechanism for the set of elementary
posteriors $\\{\pi_{k}\\}$ is chosen so as to generate an aggregated posterior
$\pi^{\star}$ that has low variance and acceptable level of bias. Moreover,
because the ultimate goal of the consensus mechanism is to mitigate the look-
ahead bias, it is imperative to choose the anchored mean $\mu$ in an
appropriate manner. As the posterior $\pi_{1}$ obtained by conditioning on the
training set $\mathcal{D}_{1}=\\{(X_{t},M_{t}):t\in T^{train}\\}$ carries the
least amount of look-ahead bias, it is reasonable to consider the discrepancy
between the mean of the aggregated posterior $\pi^{\star}$ and that of
$\pi_{1}$ as a proxy for the true bias. The optimal consensus mechanism can be
constructed from $\lambda^{\star}\in\Delta_{K}$ that solves the following
consensus optimization problem

$\begin{array}[]{cl}\min&\mathrm{Tr}(\mathrm{Cov}_{\mathcal{C}_{\lambda}(\\{\pi_{k}\\})}[\theta])\\\
\operatorname{s.t.}&\lambda\in\Delta_{K}\\\
&\ell_{\mathcal{Z}}(\mathds{E}_{\mathcal{C}_{\lambda}(\\{\pi_{k}\\})}[\theta]-\mu_{1})\leq\delta\end{array}$
(8)

for some set $\mathcal{Z}$. Notice that problem (8) shares significant
resemblance with the Markowitz’s problem: the objective function of (8)
minimizes a variance proxy, while the constraint of (8) involves the
expectation parameters.

So far, we have not yet specified how the mechanism $\mathcal{C}_{\lambda}$
depends on the parameters $\lambda$. In its most general form,
$\mathcal{C}_{\lambda}$ may depends nonlinearly on $\lambda$, and in this
case, problem (8) can be intractable. It behooves us to simplify this
dependence by imposing certain assumptions on $\Delta_{K}$ and on the
relationship between $\lambda$ and $\mathcal{C}_{\lambda}$. Towards this end,
we assume for the rest of this paper that
$\Delta_{K}=\\{\lambda\in\mathbb{R}_{+}^{K}:\sum_{k=1}^{K}\lambda_{k}=1\\}$, a
$(K-1)$-dimensional simplex. Furthermore, $\mathcal{C}_{\lambda}$ depends
explicitly on $\lambda$ through the following variational form

$\mathcal{C}_{\lambda}(\\{\pi_{k}\\})=\arg\min\limits_{\pi\in\mathcal{M}}\sum_{k\in[K]}\lambda_{k}\varphi(\pi,\pi_{k}),$
(9)

where $\varphi$ represents a measure of dissimilarity in the space of
probability measures. If $\lambda$ is a vector of $1/K$, then
$\mathcal{C}_{1/K\mathbbm{1}}(\\{\pi_{k}\\})$ is the (potentially asymmetric)
barycenter of $\\{\pi_{k}\\}$ under $\varphi$. The parameter vector $\lambda$
can now be thought of as the vector of weights corresponding to the set of
posteriors $\\{\pi_{k}\\}$, and as a consequence, the mechanism
$\mathcal{C}_{\lambda}$ outputs the weighted barycenter of $\\{\pi_{k}\\}$.

It is important to emphasize on the sharp distinction between our consensus
optimization problem (8) and the existing literature on barycenters. While the
existing literature typically focuses on studying the barycenter for a fixed
parameter $\lambda$, problem (8) searches for the best value of $\lambda$ that
meets our optimality criteria. Put differently, one can also view problem (8)
as a search problem for the best reweighting scheme of $\\{\pi_{k}\\}$ so that
their barycenter has minimal variance while being bias tolerable.

The properties of the aggregation mechanism depends fundamentally on the
choice of $\varphi$, and in what follows we explore three choices of $\varphi$
that are analytically tractable. In Section 5, we explore the mechanism using
the forward Kullback-Leibler divergence, and in Section 6, we examine the
mechanism under the type-2 Wasserstein distance. We will show that in these
two cases, the optimal consensus problem (8) is a convex optimization problem
under specific assumptions. In Section 7, we study the mechanism under the
backward Kullback-Leibler divergence, and reveal that the resulting consensus
problem is non-convex.

##  5 Optimal Forward Kullback-Leibler Consensus Mechanism

In this section, we will measure the dissimilarity between distributions using
the Kullback-Leibler (KL) divergence. The KL divergence is formally defined as
follows.

######  Definition 5.1 (Kullback-Leibler divergence).

Given two probability measures $\nu_{1}$ and $\nu_{2}$ with $\nu_{1}$ being
absolutely continuous with respect to $\nu_{2}$, the Kullback-Leibler (KL)
divergence from $\nu_{1}$ to $\nu_{2}$ is defined as
$\mathrm{KL}(\nu_{1}\\!\parallel\\!\nu_{2})\\!\triangleq\\!\mathds{E}_{\nu_{1}}\left[\log(\mathrm{d}\nu_{1}/\mathrm{d}\nu_{2})\right]$,
where $\mathrm{d}\nu_{1}/\mathrm{d}\nu_{2}$ is the Radon-Nikodym derivative of
$\nu_{1}$ with respect to $\nu_{2}$.

It is well-known that the KL divergence is non-negative, and it vanishes to
zero if and only if $\nu_{1}$ coincides with $\nu_{2}$. The KL divergence, on
the other hand, is not a distance: it even fails to be symmetric, and in
general we have
$\mathrm{KL}(\nu_{1}\parallel\nu_{2})\neq\mathrm{KL}(\nu_{2}\parallel\nu_{1})$.
The KL divergence is also a fundamental tools in statistics, information
theory (Cover and Thomas, 2012) and information geometry (Amari, 2016).

Under Assumption 3.3, the posteriors $\pi_{k}$ are Gaussian, hence, to compute
the aggregated posterior $\mathcal{C}_{\lambda}({\pi_{k}})$, we often need to
evaluate the KL divergence from, or to, a Gaussian distribution. Fortunately,
when measuring the dissimilarity between two Gaussian distributions, the KL
divergence admits a closed-form expression.

######  Lemma 5.2 (Kullback-Leibler divergence between Gaussian
distributions).

For any $(\mu_{i},\Sigma_{i})$,
$(\mu_{j},\Sigma_{j})\in\mathbb{R}^{n}\times\mathbb{S}_{++}^{n}$, we have

$\displaystyle\mathrm{KL}\big{(}\mathcal{N}_{n}(\mu_{i},\Sigma_{i})\parallel\mathcal{N}_{n}(\mu_{j},\Sigma_{j})\big{)}\triangleq(\mu_{j}-\mu_{i})^{\top}\Sigma_{j}^{-1}(\mu_{j}-\mu_{i})+\operatorname{Tr}\big{[}\Sigma_{i}\Sigma_{j}^{-1}\big{]}-\log\det(\Sigma_{i}\Sigma_{j}^{-1})-n.$

In this section we use $\varphi$ as the forward KL-divergence to aggregate the
posteriors; more specifically, for any two probability distributions $\nu_{1}$
and $\nu_{2}$, we use
$\varphi(\nu_{1},\nu_{2})=\mathrm{KL}(\nu_{1}\parallel\nu_{2})$. Note that
this choice of $\varphi$ places the posteriors $\pi_{k}$ as the second
argument in the divergence. We use $\mathcal{C}_{\lambda}^{F}$ to denote the
consensus mechanism induced by the forward KL-divergence, and
$\mathcal{C}_{\lambda}^{F}$ is formally defined as follows.

######  Definition 5.3 (Forward KL consensus mechanisms).

A mechanism $\mathcal{C}_{\lambda}^{F}$ with weight $\lambda$ is called a
forward KL consensus mechanism if for any collection of posterior
distributions $\\{\pi_{k}\\}$, we have

$\mathcal{C}_{\lambda}^{F}(\\{\pi_{k}\\})=\arg\min\limits_{\pi\in\mathcal{M}}\sum_{k\in[K]}\lambda_{k}~{}\mathrm{KL}(\pi\parallel\pi_{k}).$
(10)

The KL barycenter of the form (10) is widely used to fuse distributions and
aggregate information, and it has witnessed many successful applications in
machine learning (Claici, Yurochkin, Ghosh, and Solomon, 2020; Da, Li, Zhu,
Fan, and Fu, 2020; Li, Hlawatsch, and Djurić, 2019). The set of forward KL
mechanisms parametrized by the simplex $\Delta_{K}$ is denoted by
$\mathbb{C}^{F}\triangleq\\{\mathcal{C}_{\lambda}^{F}:\lambda\in\Delta_{K}\\}$.
For a given input $\\{\pi_{k}\\}$ and weight $\lambda\in\Delta_{K}$, computing
$\mathcal{C}_{\lambda}^{F}(\\{\pi_{k}\\})$ requires solving an infinite
dimensional optimization problem. Fortunately, when the posteriors $\pi_{k}$
are Gaussian under Assumption 3.3, we can characterize
$\mathcal{C}_{\lambda}^{F}$ analytically thanks to the following proposition.

######  Proposition 5.4 (Closed form of the consensus (Battistelli, Chisci,
Fantacci, Farina, and Graziano, 2013)).

Suppose that $\pi_{k}\sim\mathcal{N}(\mu_{k},\Sigma_{k})$ with
$\Sigma_{k}\succ 0$ for any $k\in[K]$. For any $\lambda\in\Delta_{K}$, we have
$\mathcal{C}_{\lambda}^{F}(\\{\pi_{k}\\})=\mathcal{N}(\widehat{\mu},\widehat{\Sigma})$
with

$\widehat{\mu}=(\sum_{k\in[K]}\lambda_{k}\Sigma_{k}^{-1})^{-1}\big{(}\sum_{k\in[K]}\lambda_{k}\Sigma_{k}^{-1}\mu_{k}\big{)}\in\mathbb{R}^{n},\quad\widehat{\Sigma}=(\sum_{k\in[K]}\lambda_{k}\Sigma_{k}^{-1})^{-1}\in\mathbb{S}_{++}^{n}.$
(11)

Next, we discuss the tractability of the bias constraint under the forward KL
consensus mechanism. For a given collection of posterior distributions
$\\{\pi_{k}\\}$, define the following set

$\Lambda_{\ell_{\mathcal{Z}},\delta}^{F}(\\{\pi_{k}\\})=\left\\{\lambda\in\Delta_{K}:~{}\ell_{\mathcal{Z}}(\mathds{E}_{\mathcal{C}_{\lambda}^{F}(\\{\pi_{k}\\})}[\theta]-\mu_{1})\leq\delta\right\\}$

that contains all the parameters $\lambda\in\Delta_{K}$ for which the forward
KL mechanism $\mathcal{C}_{\lambda}^{F}$ satisfies the
$(\ell_{\mathcal{Z}},\mu_{1},\delta)$-bias requirement. Unfortunately, despite
the analytical expression for the aggregated posterior in Proposition 5.4, it
remains difficult to characterize the set
$\Lambda_{\ell_{\mathcal{Z}},\delta}^{F}(\\{\pi_{k}\\})$ in its most general
form. This difficulty can be alleviated if we impose additional assumptions on
$\\{\pi_{k}\\}$ and the set $\mathcal{Z}$.

######  Proposition 5.5 ($(\ell_{\mathcal{Z}},\mu_{1},\delta)$-bias forward KL
mechanisms).

Suppose that there exists an orthogonal matrix $V\in\mathbb{R}^{n\times n}$
and a collection of diagonal matrices
$D_{k}=\operatorname{diag}(d_{k})\in\mathbb{S}_{++}^{n}$ for $k\in[K]$ such
that $\Sigma_{k}=VD_{k}V^{\top}\succ 0$ and
$\pi_{k}\sim\mathcal{N}(\mu_{k},\Sigma_{k})$ for all $k\in[K]$. If
$\mathcal{Z}=\\{z\in\mathbb{R}^{n}:\|V^{\top}z\|_{1}\leq 1\\}$, then

$\Lambda_{\ell_{\mathcal{Z}},\delta}^{F}(\\{\pi_{k}\\})=\left\\{\lambda\in\Delta_{K}:\displaystyle-\delta\leq\frac{\sum_{k\in[K]}c_{kj}\lambda_{k}}{\sum_{k\in[K]}d_{kj}^{-1}\lambda_{k}}-v_{j}^{\top}\mu_{1}\leq\delta\quad\forall
j\in[n]\right\\},$ (12)

where $d_{kj}$ is the $j$-th component of the vector $d_{k}$,
$V=[v_{1},\ldots,v_{n}]$, and

$c_{kj}=\frac{v_{j}^{\top}\mu_{k}}{d_{kj}}\quad\forall(k,j)\in[K]\times[n].$

We are now ready to state the main result of this section, which provides a
reformulation of problem (8) under conditions of Proposition 5.5.

######  Theorem 5.6 (Optimal forward KL mechanism).

Suppose that there exists an orthogonal matrix $V\in\mathbb{R}^{n\times n}$
and a collection of diagonal matrices
$D_{k}=\operatorname{diag}(d_{k})\in\mathbb{S}_{++}^{n}$ for $k\in[K]$ such
that $\Sigma_{k}=VD_{k}V^{\top}\succ 0$ and
$\pi_{k}\sim\mathcal{N}(\mu_{k},\Sigma_{k})$ for all $k\in[K]$, and suppose
that $\mathcal{Z}=\\{z\in\mathbb{R}^{n}:\|V^{\top}z\|_{1}\leq 1\\}$. Let
$\lambda^{\star}$ be the optimal solution in the variable $\lambda$ of the
following second-order cone optimization problem

$\begin{array}[]{cl}\min&\displaystyle\sum_{j\in[n]}\gamma_{j}\\\
\operatorname{s.t.}&\lambda\in\Delta_{K},~{}\gamma\in\mathbb{R}_{+}^{n}\\\
&\\!\\!\left.\begin{array}[]{l}\displaystyle\sum_{k\in[K]}c_{kj}\lambda_{k}\leq(\delta+v_{j}^{\top}\mu_{1})\sum_{k\in[K]}d_{kj}^{-1}\lambda_{k}\\\
\displaystyle\sum_{k\in[K]}c_{kj}\lambda_{k}\geq(v_{j}^{\top}\mu_{1}-\delta)\sum_{k\in[K]}d_{kj}^{-1}\lambda_{k}\\\
\left\|\begin{bmatrix}2\\\
\sum_{k\in[K]}d_{kj}^{-1}\lambda_{k}-\gamma_{j}\end{bmatrix}\right\|_{2}\leq\sum_{k\in[K]}d_{kj}^{-1}\lambda_{k}+\gamma_{j}\end{array}\right\\}\quad\forall
j\in[n],\end{array}$ (13)

where $d_{kj}$ is the $j$-th component of the vector $d_{k}$,
$V=[v_{1},\ldots,v_{n}]$, and

$c_{kj}=\frac{v_{j}^{\top}\mu_{k}}{d_{kj}}\quad\forall(k,j)\in[K]\times[n].$

Then $\mathcal{C}_{\lambda^{\star}}^{F}$ is the optimal consensus mechanism in
$\mathbb{C}^{F}$, that is, $\mathcal{C}_{\lambda^{\star}}^{F}$ is
$(\ell_{\mathcal{Z}},\mu_{1},\delta)$-bias tolerable and has minimal variance.

Problem (13) is a tractable convex optimization problem, and it can be readily
solved using off-the-shelf conic solvers such as MOSEK (MOSEK ApS, 2019). Note
that the tractability of the reformulation (13) depends crucially on the
assumption that the covariance matrices $\Sigma_{k}$ are simultaneously
diagonalizable. However, this assumption is rarely satisfied in reality. We
now discuss a pre-processing step which projects all the covariance matrices
$\Sigma_{k}$ onto a common eigenbasis. To this end, for a given orthogonal
matrix $V\in\mathbb{R}^{n\times n}$, we define the following subspace of
Gaussian distributions

$\mathcal{N}_{V}\triangleq\big{\\{}\pi\in\mathcal{M}:\exists(\mu,d)\in\mathbb{R}^{n}\times\mathbb{R}_{++}^{n},\Sigma=V\operatorname{diag}(d)V^{\top},\pi\sim\mathcal{N}(\mu,\Sigma)\big{\\}}$

that contains all nondegenerate Gaussian distributions whose covariance matrix
is diagonalizable with respect to $V$. Given any posterior
$\pi_{k}\sim\mathcal{N}(\mu_{k},\Sigma_{k})$, we add a pre-processing step to
create $\pi_{k}^{\prime}\sim\mathcal{N}(\mu_{k},\Sigma_{k}^{\prime})$ with
$\Sigma_{k}^{\prime}$ being diagonalizable with the orthogonal matrix $V$
based on the following lemma.

######  Lemma 5.7 (Projection on $\mathcal{N}_{V}$).

Let $V$ be an orthogonal matrix. Suppose that
$\pi_{k}\sim\mathcal{N}(\mu_{k},\Sigma_{k})$ for some
$\mu_{k}\in\mathbb{R}^{n}$ and $\Sigma_{k}\in\mathbb{S}_{++}^{n}$, then the
following holds:

  1. (i)

Let $d_{k}\in\mathbb{R}_{++}^{n}$ be such that
$(d_{k})_{j}=v_{j}^{\top}\Sigma_{k}v_{j}$ for $j\in[n]$. Then the distribution
$\pi_{k}^{\prime}\sim\mathcal{N}(\mu_{k},\Sigma_{k}^{\prime})$ with
$\Sigma_{k}^{\prime}=V\operatorname{diag}(d_{k})V^{\top}$ is the optimal
solution of

$\min\limits_{\pi_{k}^{\prime}\in\mathcal{N}_{V}}~{}\mathrm{KL}(\pi_{k}\parallel\pi_{k}^{\prime}).$

  2. (ii)

Let $d_{k}\in\mathbb{R}_{++}^{n}$ be such that
$(d_{k})_{j}=[V^{\top}(\Sigma_{k})^{-1}V]_{jj}$ for $j\in[n]$. Then the
distribution $\pi_{k}^{\prime}\sim\mathcal{N}(\mu_{k},\Sigma_{k}^{\prime})$
with $\Sigma_{k}^{\prime}=V\operatorname{diag}(d_{k})V^{\top}$ is the optimal
solution of

$\min\limits_{\pi_{k}^{\prime}\in\mathcal{N}_{V}}~{}\mathrm{KL}(\pi_{k}^{\prime}\parallel\pi_{k}).$

##  6 Optimal Wasserstein Consensus Mechanism

In this section, we will measure the dissimilarity between distributions using
the Wasserstein type-2 distance.

######  Definition 6.1 (Wasserstein distance).

Given two probability measures $\nu_{1}$ and $\nu_{2}$, the type-2 Wasserstein
distance between $\nu_{1}$ and $\nu_{2}$ is defined as

$\mathds{W}_{2}(\nu_{1},\nu_{2})\triangleq\min\limits_{\gamma}\left(\int_{\mathbb{R}^{n}\times\mathbb{R}^{n}}\|\theta_{1}-\theta_{2}\|_{2}^{2}\,\gamma({\rm
d}\theta_{1},{\rm d}\theta_{2})\right)^{\frac{1}{2}},$

where the minimization is taken over all couplings of $\nu_{1}$ and $\nu_{2}$.

The type-2 Wasserstein distance $\mathds{W}_{2}$ is non-negative and it
vanishes only if $\nu_{1}=\nu_{2}$. Moreover, $\mathds{W}_{2}$ is symmetric
and satisfies the triangle inequality. As a consequence, $\mathds{W}_{2}$
represents a metric on the space of probability measures (Villani, 2008, p.
94). The type-2 Wasserstein distance $\mathds{W}_{2}$ is a fundamental tool in
the theory of optimal transport, and its squared value can be interpreted as
the minimum cost of transporting the distribution $\nu_{1}$ to $\nu_{2}$, in
which the cost of moving a unit probability mass is given by the squared
Euclidean distance. The application of the Wasserstein distance has also been
extended to economics (Galichon, 2016), machine learning (Kuhn, Esfahani,
Nguyen, and Shafieezadeh-Abadeh, 2019) and statistics (Panaretos and Zemel,
2019).

As the posteriors $\pi_{k}$ are Gaussian underAssumption 3.3, we often need to
evaluate the Wasserstein distance from a Gaussian distribution. Conveniently,
the Wasserstein distance between two Gaussian distribution can be written in
an analytical form (Givens and Shortt, 1984).

######  Definition 6.2 (Wasserstein type-2 distance between Gaussian
distributions).

The Wasserstein $\mathds{W}_{2}$ distance between $\mathcal{N}(\mu,\Sigma)$
and $\mathcal{N}(\mu^{\prime},\Sigma^{\prime})$ is given by

$\mathds{W}_{2}\big{(}\mathcal{N}(\mu,\Sigma),\mathcal{N}(\mu^{\prime},\Sigma^{\prime})\big{)}=\sqrt{\|\mu-\mu^{\prime}\|_{2}^{2}+\operatorname{Tr}\big{[}\Sigma+\Sigma^{\prime}-2\big{(}(\Sigma^{\prime})^{\frac{1}{2}}\Sigma(\Sigma^{\prime})^{\frac{1}{2}}\big{)}^{\frac{1}{2}}\big{]}}.$

In this section we use $\varphi$ as the Wasserstein distance $\mathds{W}_{2}$
to aggregate the posteriors. In more details, for any two probability
distributions $\nu_{1}$ and $\nu_{2}$, we use
$\varphi(\nu_{1},\nu_{2})=\mathds{W}_{2}(\nu_{1},\nu_{2})$. We use
$\mathcal{C}_{\lambda}^{W}$ to denote the consensus mechanism induced by the
Wasserstein metric, and $\mathcal{C}_{\lambda}^{W}$ is formally defined as
follows.

######  Definition 6.3 (Wasserstein consensus mechanisms).

A mechanism $\mathcal{C}_{\lambda}^{W}$ with weight $\lambda$ is called a
Wasserstein consensus mechanism if for any collection of posterior
distributions $\\{\pi_{k}\\}\in\mathcal{N}^{K}$, we have

$\mathcal{C}_{\lambda}^{W}(\\{\pi_{k}\\})=\arg\min_{\pi\in\mathcal{N}}\sum_{k\in[K]}\lambda_{k}~{}\mathds{W}_{2}(\pi,\pi_{k})^{2}.$
(14)

Notice that by Definition 6.3, we will optimize only over the space of
Gaussian distributions in problem (14). The Wasserstein barycenter of the form
(14) is widely used to fuse distributions and aggregate information, and it
has witnessed many successful applications in machine learning, for example
distribution clustering (Ye, Wu, Wang, and Li, 2017), shape and texture
interpolation (Rabin, Peyré, Delon, and Bernot, 2012; Solomon, de Goes, Peyré,
Cuturi, Butscher, Nguyen, Du, and Guibas, 2015), and multi-target tracking
(Baum, Willett, and Hanebeck, 2015). The set of Wasserstein mechanisms
parametrized by the simplex $\Delta_{K}$ is denoted by
$\mathbb{C}^{W}\triangleq\\{\mathcal{C}_{\lambda}^{W}:\lambda\in\Delta_{K}\\}$.
When the posteriors $\pi_{k}$ are Gaussian under Assumption 3.3, we can
characterize $\mathcal{C}_{\lambda}^{W}$ semi-analytically thanks to the
following proposition.

######  Proposition 6.4 (Closed form of the consensus (Agueh and Carlier,
2011)).

Suppose that $\pi_{k}\sim\mathcal{N}(\mu_{k},\Sigma_{k})$ with
$\Sigma_{k}\succ 0$ for any $k\in[K]$. For any $\lambda\in\Delta_{K}$, we have
$\mathcal{C}_{\lambda}^{W}(\\{\pi_{k}\\})=\mathcal{N}(\widehat{\mu},\widehat{\Sigma})$
with $\widehat{\mu}=\sum_{k\in[K]}\lambda_{k}\mu_{k}\in\mathbb{R}^{n}$ and
$\widehat{\Sigma}\in\mathbb{S}_{++}^{n}$ is the unique solution of the matrix
equation

$\widehat{\Sigma}=\sum_{k\in[K]}\lambda_{k}(\widehat{\Sigma}^{\frac{1}{2}}\Sigma_{k}\widehat{\Sigma}^{\frac{1}{2}})^{\frac{1}{2}}.$
(15)

Proposition 6.4 asserts that the aggregated mean $\widehat{\mu}$ is simply the
weighted average of the posterior means $\mu_{k}$. The aggregated covariance
matrix $\widehat{\Sigma}$, unfortunately, is not available in closed-form.
Nevertheless, there exists a fast numerical routine to compute
$\widehat{\Sigma}$ for any set of input $\\{\Sigma_{k}\\}$ using fixed point
iterations (Álvarez-Esteban, del Barrio, Cuesta-Albertos, and Matrán, 2016).

In the next step, we leverage the expression of $\widehat{\mu}$ to reformulate
the bias tolerance constraint. For the collection of posterior distributions
$\\{\pi_{k}\\}$, define the following set

$\Lambda_{\ell_{\mathcal{Z}},\delta}^{W}(\\{\pi_{k}\\})=\left\\{\lambda\in\Delta_{K}:~{}\ell_{\mathcal{Z}}(\mathds{E}_{\mathcal{C}_{\lambda}^{W}(\\{\pi_{k}\\})}[\theta]-\mu_{1})\leq\delta\right\\}$

that contains the subset of $\lambda\in\Delta_{K}$ for which the Wasserstein
mechanism $\mathcal{C}_{\lambda}^{W}$ satisfies the
$(\ell_{\mathcal{Z}},\mu_{1},\delta)$-bias requirement. The next proposition
provides the expression of
$\Lambda_{\ell_{\mathcal{Z}},\delta}^{W}(\\{\pi_{k}\\})$ corresponding to the
examples of $\mathcal{Z}$ in Example 4.3.

######  Proposition 6.5 ($(\ell_{\mathcal{Z}},\mu_{1},\delta)$-bias
Wasserstein mechanisms).

The following assertions hold.

  1. (i)

If $\mathcal{Z}=\\{z\in\mathbb{R}^{n}:\|z\|_{2}\leq 1\\}$ is the Euclidean
ball, then

$\Lambda_{\ell_{\mathcal{Z}},\delta}^{W}(\\{\pi_{k}\\})=\left\\{\lambda\in\Delta_{K}:\left\|\sum_{k\in[K]}\lambda_{k}\mu_{k}-\mu_{1}\right\|_{2}\leq\delta\right\\}.$
(16)

  2. (ii)

If $\mathcal{Z}$ is a convex polyhedron
$\mathcal{Z}=\\{z\in\mathbb{R}^{n}:Az\leq b\\}$ for some matrix $A$ and vector
$b$ of appropriate dimensions, and if the linear programming

$\begin{array}[]{cl}\inf&b^{\top}w\\\
\operatorname{s.t.}&A^{\top}w=\sum_{k\in[K]}\lambda_{k}\mu_{k}-\mu_{1}\\\
&w\geq 0\end{array}$

has an optimal solution for any $\lambda\in\Delta_{K}$, then

$\Lambda_{\ell_{\mathcal{Z}},\delta}^{W}(\\{\pi_{k}\\})=\left\\{\lambda\in\Delta_{K}:\exists
w\geq
0,~{}b^{\top}w\leq\delta,~{}A^{\top}w=\left(\sum_{k\in[K]}\lambda_{k}\mu_{k}-\mu_{1}\right)\right\\}.$

  3. (iii)

If $\mathcal{Z}=\\{\mathbbm{1}_{n}\\}$, then

$\Lambda_{\ell_{\mathcal{Z}},\delta}^{W}(\\{\pi_{k}\\})=\left\\{\lambda\in\Delta_{K}:\sum_{k\in[K]}\lambda_{k}\mathbbm{1}_{n}^{\top}\mu_{k}-\mathbbm{1}_{n}^{\top}\mu_{1}\leq\delta\right\\}.$

Because the aggregated covariance matrix $\widehat{\Sigma}$ is not available
in closed form, it is intractable to devise the optimal consensus mechanism
under the Wasserstein distance in the general setting. Luckily, if we impose
the simultaneous diagonalizability assumption similar to Theorem 5.6, then the
optimal mechanism can be identified by solving a convex optimization problem,
as highlighted in the next theorem.

######  Theorem 6.6 (Optimal Wasserstein consensus mechanism).

Suppose that there exists an orthogonal matrix $V\in\mathbb{R}^{n\times n}$
and a collection of diagonal matrices
$D_{k}=\operatorname{diag}(d_{k})\in\mathbb{S}_{++}^{n}$ for $k\in[K]$ such
that $\Sigma_{k}=VD_{k}V^{\top}\succ 0$. Let $\lambda^{\star}$ be the optimal
solution of the following problem

$\begin{array}[]{cl}\min&\sum_{i,j\in[K]}\lambda_{i}\lambda_{j}\operatorname{Tr}\big{[}\Sigma_{i}^{\frac{1}{2}}\Sigma_{j}^{\frac{1}{2}}\big{]}\\\
\operatorname{s.t.}&\lambda\in\Delta_{K}\\\
&\ell_{\mathcal{Z}}\left(\sum_{k\in[K]}\lambda_{k}\mu_{k}-\mu_{1}\right)\leq\delta.\end{array}$
(17)

Then $\mathcal{C}_{\lambda^{\star}}^{W}$ is the optimal consensus mechanism in
$\mathbb{C}^{W}$, that is, $\mathcal{C}_{\lambda^{\star}}^{W}$ is
$(\ell_{\mathcal{Z}},\mu_{1},\delta)$-bias tolerable and has minimal variance.

Notice that the function

$\lambda\mapsto\sum_{i,j\in[K]}\lambda_{i}\lambda_{j}\operatorname{Tr}\big{[}\Sigma_{i}^{\frac{1}{2}}\Sigma_{j}^{\frac{1}{2}}\big{]}=\lambda^{\top}G\lambda$

is convex in $\lambda$ because the matrix
$G=(\operatorname{Tr}\big{[}\Sigma_{i}^{\frac{1}{2}}\Sigma_{j}^{\frac{1}{2}}\big{]})_{i,j\in[K]}$
is the Gram matrix of the linear reproducing kernel777More details on the
reproducing kernel Hilbert space can be found in Berlinet and Thomas-Agnan
(2004). on $\mathbb{S}_{++}^{n}$ applied to
$(\Sigma_{1}^{\frac{1}{2}},\ldots,\Sigma_{K}^{\frac{1}{2}})$, and thus $G$ is
a positive semidefinite matrix. The result of Proposition 6.5 can now be
incorporated into the reformulation (17), and lead to a convex finite-
dimensional optimization problem, which can be solved efficiently using off-
the-shelf solvers such as MOSEK (MOSEK ApS, 2019).

##  7 Optimal Backward Kullback-Leibler Consensus Mechanism

In this section, we revisit the KL divergence setting, and we consider
$\varphi$ as the backward KL-divergence, i.e., for any two probability
distributions $\nu_{1}$ and $\nu_{2}$, we use
$\varphi(\nu_{1},\nu_{2})=\mathrm{KL}(\nu_{2}\parallel\nu_{1})$. Note that
this choice of $\varphi$ places the posteriors $\pi_{k}$ as the first argument
in the divergence. We use $\mathcal{C}_{\lambda}^{B}$ to denote the consensus
mechanism induced by the backward KL-divergence, and
$\mathcal{C}_{\lambda}^{B}$ is formally defined as follows.

######  Definition 7.1 (Backward KL consensus mechanisms).

A mechanism $\mathcal{C}_{\lambda}^{B}$ with weight $\lambda$ is called a
backward KL consensus mechanism if for any collection of posterior
distributions $\\{\pi_{k}\\}$, we have

$\mathcal{C}_{\lambda}^{B}(\\{\pi_{k}\\})=\arg\min_{\pi\in\mathcal{M}}\sum_{k\in[K]}\lambda_{k}~{}\mathrm{KL}(\pi_{k}\parallel\pi).$
(18)

Given a simplex $\Delta_{K}$, we are interested in the set of backward KL
mechanisms
$\mathbb{C}^{B}\triangleq\\{\mathcal{C}_{\lambda}^{B}:\lambda\in\Delta_{K}\\}$.
Problem (18) is infinite-dimensional as the decision variable $\pi$ ranges
over the space of all probability distributions supported on $\mathbb{R}^{n}$.
Fortunately, when the posteriors $\pi_{k}$ are Gaussian under Assumption 3.3,
we can characterize $\mathcal{C}_{\lambda}^{B}$ analytically. The next
proposition states that $\mathcal{C}_{\lambda}^{B}$ is a mixture of normal
distributions.

######  Proposition 7.2 (Closed form of the consensus (Da, Li, Zhu, Fan, and
Fu, 2020, Theorem 1)).

Suppose that $\pi_{k}\sim\mathcal{N}(\mu_{k},\Sigma_{k})$ with
$\Sigma_{k}\succ 0$ for any $k\in[K]$. For any $\lambda\in\Delta_{K}$, we have
$\mathcal{C}_{\lambda}^{B}(\\{\pi_{k}\\})=\sum_{k\in[K]}\lambda_{k}\pi_{k}$, a
mixture of normal distributions.

We are now ready to state the main result of this section, which provides a
tractable reformulation of problem (8) under the backward KL consensus.

######  Theorem 7.3 (Optimal backward KL mechanism).

Suppose that $\pi_{k}\sim\mathcal{N}(\mu_{k},\Sigma_{k})$ with
$\Sigma_{k}\succ 0$ for any $k\in[K]$. Let $\lambda^{\star}$ be the optimal
solution in the variable $\lambda$ of the following problem

$\begin{array}[]{cl}\min&\sum_{k\in[K]}\lambda_{k}\operatorname{Tr}\big{[}\mu_{k}\mu_{k}^{\top}+\Sigma_{k}\big{]}-\sum_{i,j\in[K]}\lambda_{i}\lambda_{j}\mu_{i}^{\top}\mu_{j}\\\
\operatorname{s.t.}&\lambda\in\Delta_{K}\\\
&\ell_{\mathcal{Z}}\left(\sum_{k\in[K]}\lambda_{k}\mu_{k}-\mu_{1}\right)\leq\delta.\end{array}$
(19)

Then $\mathcal{C}_{\lambda^{\star}}^{B}$ is the optimal consensus mechanism in
$\mathbb{C}^{B}$, that is, $\mathcal{C}_{\lambda^{\star}}^{B}$ is
$(\ell_{\mathcal{Z}},\mu_{1},\delta)$-bias tolerable and has minimal variance.

We observe that the bias tolerance constraint in problem (19) coincide with
the bias tolerance constraint for the optimal Wasserstein consensus problem
(17). Indeed, the result of Proposition 7.2 implies that the aggregated mean
for the backward KL consensus is also a weighted average of $\mu_{k}$, and the
coincidence in the bias tolerance constraint follows. The result of
Proposition 6.5 can be reapplied in this case to reformulate the constraints
of (19) into either linear constraints or second-order cone constraints. On
the downside, notice that the objective function of (19) is a concave
quadratic function. As a consequence, solving (19) is difficult in general.
While problem (19) suffers from non-convexity, it does not require
simultaneous diagonalizability as is required in Theorem 5.6.

##  8 Simulation and Empirical Analysis

We study the look-ahead-bias and variance trade-off of our BVMI framework in
simulations and empirically. We work with a panel of asset returns data
(either simulated or empirical data). We separate the panel data into three
portions: training, testing, and an additional out-of-sample testing periods,
each having a length $T^{train},T^{test}$ and $T^{oos\text{-}test}$
respectively. For both, the simulated and empirical data, we have access to
the full panel of returns, which allows us to study various patterns of
missing data. We obtain missing data by masking entries in the panel following
certain missing mechanisms. We apply the masking procedure to the training
period, while we preserve the completeness of testing and out-of-sample
testing periods for evaluation purposes.888All optimization problems are
implemented in Python. We use an Intel Core i5 CPU (1.40 GHz) computer with
8GB of memory. Data and codes are available at
https://github.com/xuhuiustc/time-series-imputation.

Following the BVMI framework, we feed the data set with missing values to the
posterior generator module, as we have explained in Section 3.1. Recall that
we have a total of $K$ posterior distributions. The truncation times $T_{k}$’s
are set to be equispaced between $T^{train}$ and $T$ with $T_{1}=T^{train}$
and $T_{K}=T$, where $T=T^{train}+T^{test}$. We assume a non-informative flat
prior $\pi_{0}$ on the mean parameter $\theta$. The covariance matrix $\Omega$
is taken as the ground truth in the case of simulated data, or taken as the
sample covariance matrix of the complete panel in the empirical study. Our
focus is first-order problem of the uncertainty in the mean estimation.
Relative to the mean, the variance estimation is substantially more precise as
discussed among others in Kozak, Nagel, and Santosh (2020).

The generated multiple posterior distributions are used in the consensus
mechanism module to output an optimal consensus posterior, the details of
which we provide in Sections 4 to 7. More specifically, we use the following
three aggregation mechanisms:

  1. (i)

The forward KL mechanism corresponding to the reformulation (13).

  2. (ii)

The full Wasserstein mechanism corresponding to the reformulation (17),
together with the set $\mathcal{Z}$ chosen as the Euclidean ball as in
Proposition 6.5 ($i$).

  3. (iii)

A restricted Wasserstein mechanism suggested in Blanchet, Hernandez, Nguyen,
Pelger, and Zhang (2021) (Theorem 4.7), which considers the consensus
posterior as the Wasserstein barycenter of the two naive posteriors
$\pi_{1},\pi_{K}$. The use of only two posteriors gives rise to a closed-form
expression of the covariance matrix $\widehat{\Sigma}$ known as McCann’s
interpolation in McCann (1997).

We do not consider the backward KL mechanism in Section 7 as it is non-convex
and hence it is challenging to solve for large values of $K$. In order to
satisfy the simultaneous diagonalizability assumption for the forward KL
mechanism and the full Wasserstein mechanism, we project the covariance
matrices of the posteriors to the same eigenbasis. We choose the orthogonal
matrix $V$ (as in Lemma 5.7) as the eigenbasis of the covariance of the first
posterior, and we follow Lemma 5.7 (i) to project the other posteriors onto
this basis. As we have remarked before, the tolerance parameter $\delta$ in
(13) and (17) indicates how much of look-ahead bias is accepted in the
aggregated posterior $\pi^{\star}$. We select the $\delta$ parameter from the
$10$ values equi-spaced between $0$ and $\delta_{max}$, where $\delta_{max}$
is computed using (12) or (16) by setting $\lambda_{K}=1$, i.e.,
$\delta_{max}=\max_{j\in[n]}\left|c_{Kj}/d_{Kj}^{-1}-v_{j}^{\top}\mu_{1}\right|$
for the forward KL reformulation and
$\delta_{max}=\|\mu_{K}-\mu_{1}\|_{2}^{2}$ for the Wasserstein reformulation.

Corresponding to each selected $\delta$ parameter, the aggregated posterior
$\pi^{\star}$ is then output to the multiple imputation sampler module to
create multiply-imputed data sets, as we have discussed in Section 3.2. More
specifically, we follow the Conditional Expectation Bayesian Imputation
algorithm, i.e., Algorithm 2. In each round of imputation we first sample the
mean parameter $\theta$ from the aggregated posterior $\pi^{\star}$, and then
for each $t\in[T^{train}]$, we impute $Y_{t}$ by
$\theta_{Y_{t}}+\Omega_{Y_{t}X_{t}}(\Omega_{X_{t}})^{-1}(X_{t}-\theta_{X_{t}})$,
which is the conditional mean of the posterior of $Y_{t}$ given
$(\theta,\mathcal{D})$.

Our downstream task is portfolio optimization. More specifically, we construct
portfolios with the highest mean return. This captures the core of the problem
and it is straightforward to extend this to other portfolio objectives. In
more detail, the portfolio weights solve the following optimization problem on
the training period:

$\begin{array}[]{cl}\max&\mathds{E}_{\mathbb{P}^{\text{imputed}}}[w^{\top}Z]\\\
\operatorname{s.t.}&\|w\|_{2}=1,\end{array}$

where $\tilde{Z}_{t}$ denotes the imputed value of $Z_{t}$ and
$\mathbb{P}^{\text{imputed}}=\sum_{t\in[T^{train}]}\delta_{\tilde{Z}_{t}}/T^{train}$
is the empirical distribution supported on the imputed data. The analytical
solution equals

$w=\displaystyle\frac{\frac{1}{T^{train}}\sum_{t\in[T^{train}]}\tilde{Z}_{t}}{\|\frac{1}{T^{train}}\sum_{t\in[T^{train}]}\tilde{Z}_{t}\|_{2}}.$

The estimated portfolio weights are used on the testing period, which results
in the average portfolio return of

$R_{test}=\frac{1}{T^{test}}\sum_{t\in[T^{train}+1,T]}w^{\top}Z_{t}.$

Similarly, when the portfolio weights $w$ are applied on the out-of-sample
testing period, we obtain the averaged portfolio return

$R_{oos\text{-}test}=\frac{1}{T^{oos\text{-}test}}\sum_{t\in[T+1,T+T^{oos\text{-}test}]}w^{\top}Z_{t}.$

The difference in the averaged portfolio return between the two periods
defines our regret measure

$\Delta R=R_{test}-R_{oos\text{-}test}.$

Recall that $\mathcal{D}$ is the realization of observed data and missing
masks. Conditional on $\mathcal{D}$, $\Delta R$ is a random variable where the
randomness comes from multiple imputations. Thus we have the conditional
moments $\mathbb{E}[\Delta R|\mathcal{D}]$ and $\mathbb{V}\mathrm{ar}[\Delta
R|\mathcal{D}]$, whose randomness can in turn be averaged out over
realizations of $\mathcal{D}$. We are interested in finding a $\delta$ that
minimizes the “expected conditional mean squared error”, or ECMSE

$\max\\{\mathbb{E}[\Delta R],0\\}^{2}+\mathbb{E}[\mathbb{V}\mathrm{ar}[\Delta
R|\mathcal{D}]].$

Observe that this ECMSE shares the same bias-variance trade off structure as
problem (8) (with $\mathcal{Z}$ chosen according to Example 4.3 ($i$)) with
the exception that we consider a one-sided bias in the above formulation as an
investor’s objective is to avoid a negative regret or disappointment. We
denote the first term above as ECBias2 and the second term above as ECVar. In
the following we study the three measure ECMSE, ECBias2 and ECVar as $\delta$
varies, under four different missing mechanisms that includes missing at
random, a dependent block missing structure, and when the missing pattern
depends on realizations of returns (noticeably the last missing pattern does
not satisfy the ignorability assumption, i.e., our Assumption 3.2).

###  8.1 Simulation

In our simulation we generate returns for $n=10$ stocks based on a factor
model and study the effect of data imputation for different missing patterns.
We consider the time horizons $T^{train}=T^{test}=100$,
$T^{oos\text{-}test}=1000$. The large size of testing set is designed to
amplify the look-ahead-bias. We generate the returns vector from a factor
model, namely, the vectors $Z_{t}$,
$t=1,\ldots,T^{train}+T^{test}+T^{oos\text{-}test}$ are generated i.i.d. from
$\mathcal{N}_{n}(\theta,\Omega)$, where for each $i=1,\ldots,n$,
$\theta_{i}=0.2\beta_{i}+\alpha_{i}$, $\beta_{i}=1$, $\alpha_{i}$ is equi-
spaced between $-0.3$ and $0.3$, and $\Omega=\beta\beta^{\top}+I_{n}$.

We consider four types of missing mechanisms.

  * •

Missing completely at random – we mask each entry in the training period as
missing independently at random with a fixed probability $50\%$.

  * •

Missing at random – for each stock $i$, we sample a Bernoulli variable $S_{i}$
with success probability $50\%$. If $S_{i}=1$, we mask the entries
corresponding to stock $i$ as missing independently at random with probability
$50\%$, otherwise we mask the entries corresponding to stock $i$ as missing
independently at random with probability $70\%$.

  * •

Block missing – we mask the first $30\%$ of the training period as missing.

  * •

Missing by value – we mask an entry in the training period as missing if it
has an absolute value greater than $0.3$.

Figure 3: Model comparison on simulated data

(a) Missing completely at random

(b) Missing at random

(c) Block missing

(d) Missing by value

These figures compares the expected conditional mean squared error (ECMSE) for
different imputation schemes on simulated data. We impute the missing data
using the forward KL mechanism (blue solid line with circle), full Wasserstein
mechanism (green dash-dot line with square) and restricted Wasserstein
mechanism (dashed orange line with pentagon) for different weights $\delta$,
which control the trade-off between look-ahead-bias and variance.

For each missing mechanism, we run $500$ independent simulations, where in
each simulated sample the imputation is repeated for $10$ times, and we carry
out the downstream task on each imputed data set. Our main focus is the effect
of the parameter $\delta$ for the three aggregation mechanisms, forward KL
mechanism, full Wasserstein mechanism and restricted Wasserstein mechanism.
Figure 3 shows the comparison of the ECMSE for the three mechanisms and the
four missing patterns. Importantly, for all missing patterns and mechanisms
the expected conditional mean squared error is minimized for a non-trivial
value of $\delta$. This means the special cases of no-look-ahead bias
($\delta=0$) or using all the data with maximum look-ahead-bias
($\delta=\delta_{max}$) are not optimal. Moreover, either the forward KL
mechanism or the full Wasserstein mechanism attains a smaller minimum ECMSE
than the restricted Wasserstein mechansim, showing the power of using more
layers of information in the aggregation mechanism.

Figure 4: Look-ahead-bias and variance with forward KL mechanism on simulated
data

(a) Missing completely at random

(b) Missing at random

(c) Block missing

(d) Missing by value

These figures show the look-ahead-bias measured by the regret ECBias2, the
variance ECVar and the expected conditional mean squared error ECMSE for
simulated data with the forward KL mechanism for different weights $\delta$,
which control the trade-off between look-ahead-bias and variance.

Figure 4 sheds further light on the underlying trade-off and shows the ECBias2
and ECVar for the forward KL mechanism. In all cases the look-ahead-bias
measured by ECBias2 is increasing in $\delta$, while the variance measured by
ECVar decreases. This is the reason why ECMSE can achieve its optimal value
for non-trivial values of $\delta$. Figures A.1 and A.2 in the Appendix show
the corresponding results for the full and restricted Wasserstein mechanism
with the same monotonicity results.

###  8.2 Empirical Application

In our empirical application we consider two different panels of daily stock
returns that cover different time periods. Our first data set is a
representative sample of $n=10$ large capitalization stocks999The stocks are
‘DVA’, ‘DE’, ‘DAL’,‘ESRX’, ‘XOM’, ‘FFIV’, ‘FB’, ‘FAST’, ‘ESS’ and ‘EL’. from
01/01/2015 to 01/01/2017. Our second panel consists of the $n=10$ standard
industry portfolios from Kenneth French’s data library from 01/03/2017 to
01/02/2019.101010The industry portfolios are ‘NoDur’, ‘Durbl’, ‘Manuf’,
‘Enrgy’, ‘HiTec’, ‘Telcm’, ‘Shops’, ‘Hlth’, ‘Utils’ and ‘Other’, downloaded
from
https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html. We
consider $T^{train}=200$, $T^{test}=100$ and $T^{oos\text{-}test}=100$. All
returns are annualized. Similar to before, we consider four types of missing
mechanisms:

Figure 5: Model comparison for individual stocks

(a) Missing completely at random

(b) Missing at random

(c) Block missing

(d) Missing by value

These figures compares the expected conditional mean squared error (ECMSE) for
different imputation schemes for individual stock returns. We impute the
missing data using the forward KL mechanism (blue solid line with circle),
full Wasserstein mechanism (green dash-dot line with square) and restricted
Wasserstein mechanism (dashed orange line with pentagon) for different weights
$\delta$, which control the trade-off between look-ahead-bias and variance.
The sample are daily returns of 10 representative large capitalization stocks
from 01/01/2015 to 01/01/2017. Figure 6: Look-ahead-bias and variance with
forward KL mechanism for individual stocks

(a) Missing completely at random

(b) Missing at random

(c) Block missing

(d) Missing by value

These figures show the look-ahead-bias measured by the regret ECBias2, the
variance ECVar and the expected conditional mean squared error ECMSE for
individual stocks with the forward KL mechanism for different weights
$\delta$, which control the trade-off between look-ahead-bias and variance.
The sample are daily returns of 10 representative large capitalization stocks
from 01/01/2015 to 01/01/2017. Figure 7: Model comparison for industry
portfolios

(a) Missing completely at random

(b) Missing at random

(c) Block missing

(d) Missing by value

These figures compares the expected conditional mean squared error (ECMSE) for
different imputation schemes for portfolio returns. We impute the missing data
using the forward KL mechanism (blue solid line with circle), full Wasserstein
mechanism (green dash-dot line with square) and restricted Wasserstein
mechanism (dashed orange line with pentagon) for different weights $\delta$,
which control the trade-off between look-ahead-bias and variance. The sample
are daily returns of 10 industry portfolios from 01/03/2017 to 01/02/2019.
Figure 8: Look-ahead-bias and variance with forward KL mechanism for industry
portfolios

(a) Missing completely at random

(b) Missing at random

(c) Block missing

(d) Missing by value

These figures show the look-ahead-bias measured by the regret ECBias2, the
variance ECVar and the expected conditional mean squared error ECMSE for
industry portfolios with the forward KL mechanism for different weights
$\delta$, which control the trade-off between look-ahead-bias and variance.
The sample are daily returns of 10 industry portfolios from 01/03/2017 to
01/02/2019.

  * •

Missing completely at random – we mask each entry in the training period as
missing independently at random with a fixed probability $40\%$. The data is
fixed as the first $T^{train}+T^{test}+T^{oos\text{-}test}=400$ dates from the
original data set. We repeat $50$ independent simulations, where in each
sample we first generate the missing masks and then multiply-impute $100$
times.

  * •

Missing at random – for each industry portfolio $i$, we sample a Bernoulli
variable $S_{i}$ with success probability $50\%$. If $S_{i}=1$, we mask the
entries corresponding to $i$ as missing independently at random with
probability $20\%$, otherwise we mask the entries corresponding to $i$ as
missing independently at random with probability $50\%$. The rest of the
procedure is the same as the previous case.

  * •

Block missing – we mask the first $40\%$ of the training period as missing. We
conduct $50$ independent simulations. In the $k$-th sample, we use the
$T^{train}+T^{test}+T^{oos\text{-}test}=400$ days of data starting from the
$k$-th date in the original data set, and we multiply-impute $100$ times.

  * •

Missing by value – we mask an entry in the training period as missing if its
absolute value is greater than $0.03$. The rest of the procedure is the same
as the previous case.

Figure 5 compares the ECMSE as function of $\delta$ for the three aggregation
mechanisms (forward KL mechanism, full Wasserstein mechanism and restricted
Wasserstein mechanism) for individual stock returns. First, we observe that in
most cases a non-trivial $\delta$ is optimal, that is, neither avoiding any
look-ahead-bias ($\delta=0$) nor using all the data with full look-ahead-bias
($\delta=\delta_{max}$) minimizes the ECMSE. Second, the optimal value of
$\delta$ depends on the missing pattern. Third, while the three mechanism
provide somewhat similar results and trade-offs, the more flexible forward KL
and full Wasserstein mechanism dominate the restricted Wasserstein mechanism.
Figure 6 provides further details for the trade-off between ECBsias2 and ECVar
for the forward KL mechanism. The results for the Wasserstein mechanism are
similar and collected in Figures A.3 and A.4 in the Appendix. The monotonic
dependency on $\delta$ is as expected. In most cases the variance is only
marginally decreased by using more data, while the increase in the bias seems
to drive the trade-off.

Figures 7, 8, A.5 and A.6 show the corresponding results for the industry
portfolios. While the general findings are the same as for the individual
stocks, they highlight the important dependency on the missing pattern. In the
case of missing at random a small value of $\delta$ is optimal as the look-
ahead-bias is increasing much faster than the variance. However, for more
complex missing patterns, where the returns are missing in blocks or for large
realizations, the decrease in variance becomes more dominant and hence leads
to a larger optimal value of $\delta$.

##  9 Conclusion

In this paper we study the implications of imputation in a downstream out-of-
sample task for financial time-series data. We propose the Bayesian framework
to introduce and precisely formulate the notion of look-ahead-bias involved in
imputation procedures. We provide a solution for the optimal trade-off between
look-ahead-bias and variance by a consensus posterior that fuses a number of
individual posteriors corresponding to different layers of future information,
and we use the (forward and backward) KL divergence and Wasserstein distance
to derive tractable reformulations for finding the optimal consensus
posterior. Our theoretical results focus on the fundamental case of estimating
mean returns under simplifying assumptions in order to clearly illustrate the
novel conceptual ideas. However, the stylized features we adopt have a solid
motivation in finance and can be easily relaxed. For example, an extension of
our analysis to residuals following an AR(p) process is straightforward
(Little and Rubin (2002)).

Last but not least, we want to emphasize that the impact of look-ahead-bias is
important in financial applications and even relatively small improvements are
valuable because their impact can be largely magnified in the downstream task.
An investment strategy that applies leverage or derivatives can substantially
increase the numbers reported in this paper. We are the first to identify and
study the problem of look-ahead-bias in the context of data imputation. As the
use of increasingly large data sets is consistently gaining relevance in
finance, while most of these data sets have missing values, it is crucial to
understand the implications of how to deal with missing values. We are
confident that we have “opened a floodgate” and that more follow-up papers
will focus on the issue of look-ahead-bias quantification, which we believe is
interesting, relevant and largely overlooked.

## References

  * (1)
  * Agueh and Carlier (2011) Agueh, M., and G. Carlier (2011): “Barycenters in the Wasserstein Space,” _SIAM Journal on Mathematical Analysis_ , 43(2), 904–924.
  * Álvarez-Esteban, del Barrio, Cuesta-Albertos, and Matrán (2016) Álvarez-Esteban, P. C., E. del Barrio, J. Cuesta-Albertos, and C. Matrán (2016): “A fixed-point approach to barycenters in Wasserstein space,” _Journal of Mathematical Analysis and Applications_ , 441(2), 744–762.
  * Amari (2016) Amari, S. (2016): _Information Geometry and Its Applications_. Springer.
  * Audigier, Husson, and Josse (2016) Audigier, V., F. Husson, and J. Josse (2016): “A principal component method to impute missing values for mixed data,” _Advances in Data Analysis and Classification_ , 10(1), 5–26.
  * Bai and Ng (2020) Bai, J., and S. Ng (2020): “Matrix Completion, Counterfactuals, and Factor Analysis of Missing Data,” arXiv:1910.06677.
  * Battistelli, Chisci, Fantacci, Farina, and Graziano (2013) Battistelli, G., L. Chisci, C. Fantacci, A. Farina, and A. Graziano (2013): “Consensus CPHD Filter for Distributed Multitarget Tracking,” _IEEE Journal of Selected Topics in Signal Processing_ , 7(3), 508–520.
  * Baum, Willett, and Hanebeck (2015) Baum, M., P. K. Willett, and U. D. Hanebeck (2015): “On Wasserstein Barycenters and MMOSPA Estimation,” _IEEE Signal Processing Letters_ , 22(10), 1511–1515.
  * Berlinet and Thomas-Agnan (2004) Berlinet, A., and C. Thomas-Agnan (2004): _Reproducing Kernel Hilbert Space in Probability and Statistics_. Kluwer Academic Publishers.
  * Blanchet, Hernandez, Nguyen, Pelger, and Zhang (2021) Blanchet, J., F. Hernandez, V. A. Nguyen, M. Pelger, and X. Zhang (2021): “Time-Series Imputation with Wasserstein Interpolation for Optimal Look-Ahead-Bias and Variance Tradeoff,” arXiv:2102.12736.
  * Bryzgalova, Lettau, Lerner, and Pelger (2022) Bryzgalova, S., M. Lettau, S. Lerner, and M. Pelger (2022): “Asset Pricing with Missing Data,” _Working paper_.
  * Buchholz, Ahfock, and Richardson (2019) Buchholz, A., D. Ahfock, and S. Richardson (2019): “Distributed Computation for Marginal Likelihood based Model Choice,” arXiv:1910.04672.
  * Cai, Candès, and Shen (2010) Cai, J.-F., E. J. Candès, and Z. Shen (2010): “A Singular Value Thresholding Algorithm for Matrix Completion,” _SIAM Journal on Optimization_ , 20(4), 1956–1982.
  * Cambanis, Huang, and Simons (1981) Cambanis, S., S. Huang, and G. Simons (1981): “On the theory of elliptically contoured distributions,” _Journal of Multivariate Analysis_ , 11(3), 368–385.
  * Candès and Recht (2009) Candès, E. J., and B. Recht (2009): “Exact Matrix Completion via Convex Optimization,” _Foundations of Computational Mathematics_ , 9(6), 717.
  * Candès and Tao (2010) Candès, E. J., and T. Tao (2010): “The Power of Convex Relaxation: Near-Optimal Matrix Completion,” _IEEE Transactions on Information Theory_ , 56(5), 2053–2080.
  * Cao, Wang, Li, Zhou, Li, and Li (2018) Cao, W., D. Wang, J. Li, H. Zhou, Y. Li, and L. Li (2018): “BRITS: Bidirectional Recurrent Imputation for Time Series,” in _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , pp. 6776–6786.
  * Che, Purushotham, Cho, Sontag, and Liu (2018) Che, Z., S. Purushotham, K. Cho, D. Sontag, and Y. Liu (2018): “Recurrent Neural Networks for Multivariate Time Series with Missing Values,” _Nature Scientific Reports_ , 8(1), 6085.
  * Chen, Pelger, and Zhu (2019) Chen, L., M. Pelger, and J. Zhu (2019): “Deep Learning in Asset Pricing,” _Working paper_.
  * Chen, Fan, Ma, and Yan (2019) Chen, Y., J. Fan, C. Ma, and Y. Yan (2019): “Inference and uncertainty quantification for noisy matrix completion,” _Proceedings of the National Academy of Sciences_ , 116(46), 22931–22937.
  * Choi, Bahadori, Schuetz, Stewart, and Sun (2016) Choi, E., M. T. Bahadori, A. Schuetz, W. F. Stewart, and J. Sun (2016): “Doctor AI: Predicting Clinical Events via Recurrent Neural Networks,” vol. 56 of _Proceedings of Machine Learning Research_ , pp. 301–318.
  * Claici, Yurochkin, Ghosh, and Solomon (2020) Claici, S., M. Yurochkin, S. Ghosh, and J. Solomon (2020): “Model Fusion with Kullback-Leibler Divergence,” in _Proceedings of the 37th International Conference on Machine Learning_ , pp. 2038–2047.
  * Clogg, Rubin, Schenker, Schultz, and Weidman (1991) Clogg, C. C., D. B. Rubin, N. Schenker, B. Schultz, and L. Weidman (1991): “Multiple Imputation of Industry and Occupation Codes in Census Public-Use Samples Using Bayesian Logistic Regression,” _Journal of the American Statistical Association_ , 86(413), 68–78.
  * Cover and Thomas (2012) Cover, T., and J. Thomas (2012): _Elements of Information Theory_. John Wiley & Sons.
  * Da, Li, Zhu, Fan, and Fu (2020) Da, K., T. Li, Y. Zhu, H. Fan, and Q. Fu (2020): “Kullback-Leibler Averaging for Multitarget Density Fusion,” in _16th International Conference on Distributed Computing and Artificial Intelligence_ , pp. 253–261.
  * Dai, Pollock, and Roberts (2021) Dai, H., M. Pollock, and G. Roberts (2021): “Bayesian Fusion: Scalable unification of distributed statistical analyses,” arXiv:2102.02123.
  * de Leeuw and van der Heijden (1988) de Leeuw, J., and P. G. M. van der Heijden (1988): “Correspondence analysis of incomplete contingency tables,” _Psychometrika_ , 53(2), 223–233.
  * DeMiguel, Garlappi, and Uppal (2007) DeMiguel, V., L. Garlappi, and R. Uppal (2007): “Optimal Versus Naive Diversification: How Inefficient is the $1/N$ Portfolio Strategy?,” _The Review of Financial Studies_ , 22(5), 1915–1953.
  * Dempster, Laird, and Rubin (1977) Dempster, A. P., N. M. Laird, and D. B. Rubin (1977): “Maximum Likelihood from Incomplete Data via the EM Algorithm,” _Journal of the Royal Statistical Society. Series B_ , 39(1), 1–38.
  * Fang, Kotz, and Ng (1990) Fang, K.-T., S. Kotz, and K. W. Ng (1990): _Symmetric Multivariate and Related Distributions_. Chapman & Hall.
  * Friedman, Huang, Zhang, and Cao (2012) Friedman, C., J. Huang, Y. Zhang, and W. Cao (2012): “Missing value imputation via copula and transformation methods, with applications to financial and economic data,” _International Journal of Data Analysis Techniques and Strategies_ , 4(4), 315–339.
  * Galichon (2016) Galichon, A. (2016): _Optimal Transport Methods in Economics_. Princeton University Press.
  * Giglio, Liao, and Xiu (2021) Giglio, S., Y. Liao, and D. Xiu (2021): “Thousands of Alpha Tests,” _Review of Financial Studies_ , 34(7), 3456–3496.
  * Givens and Shortt (1984) Givens, C., and R. Shortt (1984): “A class of Wasserstein metrics for probability distributions,” _The Michigan Mathematical Journal_ , 31(2), 231–240.
  * Gu, Kelly, and Xiu (2020) Gu, S., B. T. Kelly, and D. Xiu (2020): “Empirical Asset Pricing via Machine Learning,” _Review of Financial Studies_ , 33(5), 2223–2273.
  * Halpin (2016) Halpin, B. (2016): “Multiple Imputation for Categorical Time Series,” _The Stata Journal_ , 16(3), 590–612.
  * Honaker and King (2010) Honaker, J., and G. King (2010): “What to do About Missing Values in Time Series Cross-Section Data,” _American Journal of Political Science_ , 54(3), 561–581.
  * Hopke, Liu, and Rubin (2001) Hopke, P. K., C. Liu, and D. B. Rubin (2001): “Multiple Imputation for Multivariate Data with Missing and Below-Threshold Measurements: Time-Series Concentrations of Pollutants in the Arctic,” _Biometrics_ , 57(1), 22–33.
  * Keshavan, Montanari, and Oh (2010) Keshavan, R. H., A. Montanari, and S. Oh (2010): “Matrix Completion From a Few Entries,” _IEEE Transactions on Information Theory_ , 56(6), 2980–2998.
  * Kozak, Nagel, and Santosh (2020) Kozak, S., S. Nagel, and S. Santosh (2020): “Shrinking the Cross-Section,” _Journal of Financial Economics_ , 135(2), 271–292.
  * Kreindler and Lumsden (2006) Kreindler, D. M., and C. J. Lumsden (2006): “The effects of the irregular sample and missing data in time series analysis,” _Nonlinear Dynamics, Psychology, and Life Sciences_ , 10(2), 187–214.
  * Kuhn, Esfahani, Nguyen, and Shafieezadeh-Abadeh (2019) Kuhn, D., P. M. Esfahani, V. A. Nguyen, and S. Shafieezadeh-Abadeh (2019): “Wasserstein Distributionally Robust Optimization: Theory and Applications in Machine Learning,” _INFORMS TutORials in Operations Research_ , pp. 130–166.
  * Lan, Xu, Ma, and Li (2020) Lan, Q., X. Xu, H. Ma, and G. Li (2020): “Multivariable data imputation for the analysis of incomplete credit data,” _Expert Systems with Applications_ , 141, 112926.
  * Li, Hlawatsch, and Djurić (2019) Li, T., F. Hlawatsch, and P. M. Djurić (2019): “Cardinality-Consensus-Based PHD Filtering for Distributed Multitarget Tracking,” _IEEE Signal Processing Letters_ , 26(1), 49–53.
  * Lipton, Kale, and Wetzel (2016) Lipton, Z. C., D. Kale, and R. Wetzel (2016): “Directly Modeling Missing Data in Sequences with RNNs: Improved Classification of Clinical Time Series,” vol. 56 of _Proceedings of Machine Learning Research_ , pp. 253–270.
  * Little and Rubin (2002) Little, R., and D. Rubin (2002): _Statistical Analysis with Missing Data_. Wiley.
  * Lobo, Vandenberghe, Boyd, and Lebret (1998) Lobo, M. S., L. Vandenberghe, S. Boyd, and H. Lebret (1998): “Applications of second-order cone programming,” _Linear Algebra and its Applications_ , 284(1), 193 – 228.
  * Mayer, Josse, Tierney, and Vialaneix (2019) Mayer, I., J. Josse, N. Tierney, and N. Vialaneix (2019): “R-miss-tastic: a unified platform for missing values methods and workflows,” arXiv:1908.04822.
  * Mazumder, Hastie, and Tibshirani (2010) Mazumder, R., T. Hastie, and R. Tibshirani (2010): “Spectral Regularization Algorithms for Learning Large Incomplete Matrices,” _Journal of Machine Learning Research_ , 11(80), 2287–2322.
  * McCann (1997) McCann, R. J. (1997): “A Convexity Principle for Interacting Gases,” _Advances in Mathematics_ , 128(1), 153–179.
  * MOSEK ApS (2019) MOSEK ApS (2019): _MOSEK Optimizer API for Python 9.2.10_.
  * Panaretos and Zemel (2019) Panaretos, V. M., and Y. Zemel (2019): “Statistical Aspects of Wasserstein Distances,” _Annual Review of Statistics and Its Application_ , 6(1), 405–431.
  * Rabin, Peyré, Delon, and Bernot (2012) Rabin, J., G. Peyré, J. Delon, and M. Bernot (2012): “Wasserstein Barycenter and Its Application to Texture Mixing,” in _Scale Space and Variational Methods in Computer Vision_ , ed. by A. M. Bruckstein, B. M. ter Haar Romeny, A. M. Bronstein, and M. M. Bronstein, pp. 435–446.
  * Rubin (1987) Rubin, D. B. (1987): _Multiple Imputation for Nonresponse in Surveys_. Wiley.
  * Schervish (1995) Schervish, M. J. (1995): _Theory of Statistics_. Springer.
  * Scott, Blocker, Bonassi, Chipman, George, and McCulloch (2016) Scott, S. L., A. W. Blocker, F. V. Bonassi, H. A. Chipman, E. I. George, and R. E. McCulloch (2016): “Bayes and Big Data: The Consensus Monte Carlo Algorithm,” _International Journal of Management Science and Engineering Management_ , 11, 78–88.
  * Solomon, de Goes, Peyré, Cuturi, Butscher, Nguyen, Du, and Guibas (2015) Solomon, J., F. de Goes, G. Peyré, M. Cuturi, A. Butscher, A. Nguyen, T. Du, and L. Guibas (2015): “Convolutional Wasserstein Distances: Efficient Optimal Transportation on Geometric Domains,” _ACM Trans. Graph._ , 34(4).
  * Tsay (2010) Tsay, R. (2010): _Analysis of Financial Time Series_. Wiley.
  * van Buuren (2012) van Buuren, S. (2012): _Flexible Imputation of Missing Data_. CRC Press.
  * Villani (2008) Villani, C. (2008): _Optimal Transport: Old and New_. Springer Science & Business Media.
  * Xiong and Pelger (2019) Xiong, R., and M. Pelger (2019): “Large Dimensional Latent Factor Modeling with Missing Observations and Applications to Causal Inference,” arXiv:1910.08273.
  * Ye, Wu, Wang, and Li (2017) Ye, J., P. Wu, J. Z. Wang, and J. Li (2017): “Fast Discrete Distribution Clustering Using Wasserstein Barycenter With Sparse Support,” _Transactions on Signal Processing_ , 65(9), 2317–2332.
  * Yoon, Jordon, and van der Schaar (2018) Yoon, J., J. Jordon, and M. van der Schaar (2018): “GAIN: Missing Data Imputation using Generative Adversarial Nets,” vol. 80 of _Proceedings of Machine Learning Research_ , pp. 5689–5698.

##  Appendix A Proof of Results

###### Proof of Lemma 3.4.

Let $\xi$ by any non-zero vector of dimension $n$. From Assumption 3.1, there
exists $t\in[T_{k}]$, such that $\xi_{X_{t}}=\mathcal{P}_{M_{t}}^{\perp}(\xi)$
is non-zero. Since $\Omega_{X_{t}}$ is positive definite, we have
$\xi_{X_{t}}^{\top}\Omega_{X_{t}}^{-1}\xi_{X_{t}}>0$. Therefore

$\xi^{\top}\left((\mathcal{P}_{M_{t}}^{\perp})^{-1}(\Omega_{X_{t}}^{-1})\right)\xi>0.$

Thus

$\xi^{\top}\left(\sum_{t\in[T_{k}]}(\mathcal{P}_{M_{t}}^{\perp})^{-1}(\Omega_{X_{t}}^{-1})\right)\xi>0.$

Since $\xi$ is arbitrary, we have
$\sum_{t\in[T_{k}]}(\mathcal{P}_{M_{t}}^{\perp})^{-1}(\Omega_{X_{t}}^{-1})$ is
positive definite. Thus $\Sigma_{k}$ is positive definite. ∎

###### Proof of Proposition 5.5.

From the expression of the aggregated mean $\widehat{\mu}$ in Proposition 5.4,
we have

$\displaystyle\widehat{\mu}$
$\displaystyle=\big{(}\sum_{k\in[K]}\lambda_{k}\Sigma_{k}^{-1}\big{)}^{-1}\big{(}\sum_{k\in[K]}\lambda_{k}\Sigma_{k}^{-1}\mu_{k}\big{)}$
$\displaystyle=\big{(}\sum_{k\in[K]}\lambda_{k}VD_{k}^{-1}V^{\top}\big{)}^{-1}\big{(}\sum_{k\in[K]}\lambda_{k}VD_{k}^{-1}V^{\top}\mu_{k}\big{)}$
$\displaystyle=V\big{(}\sum_{k\in[K]}\lambda_{k}D_{k}^{-1}\big{)}^{-1}V^{\top}V\big{(}\sum_{k\in[K]}\lambda_{k}D_{k}^{-1}V^{\top}\mu_{k}\big{)}$
$\displaystyle=V\big{(}\sum_{k\in[K]}\lambda_{k}D_{k}^{-1}\big{)}^{-1}\big{(}\sum_{k\in[K]}\lambda_{k}D_{k}^{-1}V^{\top}\mu_{k}\big{)},$

where the second equality follows from the eigendecomposition of $\Sigma_{k}$
and the last equality follows because $V$ is an orthogonal matrix and thus
$V^{\top}V=I_{n}$. By exploiting the definition of the set $\mathcal{Z}$, we
have

$\ell_{\mathcal{Z}}(\mathds{E}_{\mathcal{C}_{\lambda}^{F}(\\{\pi_{k}\\})}[\theta]-\mu_{1})=\ell_{\mathcal{Z}}(\widehat{\mu}-\mu_{1})=\sup\limits_{z\in\mathcal{Z}}z^{\top}(\widehat{\mu}-\mu_{1})=\sup\limits_{z:\|V^{\top}z\|_{1}\leq
1}z^{\top}VV^{\top}(\widehat{\mu}-\mu_{1})=\|V^{\top}(\widehat{\mu}-\mu_{1})\|_{\infty},$

where the third equality follows from the fact that $VV^{\top}=I_{n}$, and the
last equality follows by the properties of the dual norm pairs between
$\|\cdot\|_{1}$ and $\|\cdot\|_{\infty}$. Combine with the expression for
$\widehat{\mu}$ derived previously, we have

$\displaystyle\ell_{\mathcal{Z}}(\widehat{\mu}-\mu_{1})$
$\displaystyle=\|\big{(}\sum_{k\in[K]}\lambda_{k}D_{k}^{-1}\big{)}^{-1}\big{(}\sum_{k\in[K]}\lambda_{k}D_{k}^{-1}V^{\top}\mu_{k}\big{)}-V^{\top}\mu_{1}\|_{\infty}$
$\displaystyle=\max_{j\in[n]}\left|\frac{\sum_{k\in[K]}c_{kj}\lambda_{k}}{\sum_{k\in[K]}d_{kj}^{-1}\lambda_{k}}-v_{j}^{\top}\mu_{1}\right|.$

Therefore the following holds

$\\{\lambda\in\Delta_{K}:\ell_{\mathcal{Z}}(\mathds{E}_{\mathcal{C}_{\lambda}^{F}(\\{\pi_{k}\\})}[\theta]-\mu_{1})\leq\delta\\}=\left\\{\lambda\in\Delta_{K}:\displaystyle-\delta\leq\frac{\sum_{k\in[K]}c_{kj}\lambda_{k}}{\sum_{k\in[K]}d_{kj}^{-1}\lambda_{k}}-v_{j}^{\top}\mu_{1}\leq\delta\quad\forall
j\in[n]\right\\},$

which is a feasible set prescribed using one simplex constraint and a
collection of fractional linear constraints. This completes the proof. ∎

###### Proof of Theorem 5.6.

By exploiting Proposition 5.4, we have
$\mathcal{C}_{\lambda}^{F}(\\{\pi_{k}\\})=\mathcal{N}(\widehat{\mu},\widehat{\Sigma})$,
therefore problem (8) is equivalent to

$\begin{array}[]{cl}\min&\operatorname{Tr}\big{[}\widehat{\Sigma}\big{]}\\\
\operatorname{s.t.}&\lambda\in\Delta_{K}\\\
&\widehat{\Sigma}=\big{(}\sum_{k\in[K]}\lambda_{k}\Sigma_{k}^{-1}\big{)}^{-1}\in\mathbb{S}_{++}^{n},~{}\widehat{\mu}=\widehat{\Sigma}\big{(}\sum_{k\in[K]}\lambda_{k}\Sigma_{k}^{-1}\mu_{k}\big{)}\in\mathbb{R}^{n}\\\
&\ell_{\mathcal{Z}}(\widehat{\mu}-\mu_{1})\leq\delta.\end{array}$ (A.1)

By exploiting the diagonalization of $\Sigma_{k}$, we have

$\displaystyle\operatorname{Tr}\big{[}\widehat{\Sigma}\big{]}$
$\displaystyle=\operatorname{Tr}\big{[}\big{(}\sum_{k\in[K]}\lambda_{k}\Sigma_{k}^{-1}\big{)}^{-1}\big{]}=\operatorname{Tr}\big{[}\big{(}\sum_{k\in[K]}\lambda_{k}VD_{k}^{-1}V^{\top}\big{)}^{-1}\big{]}$
$\displaystyle=\operatorname{Tr}\big{[}V(\sum_{k\in[K]}\lambda_{k}D_{k}^{-1})^{-1}V^{\top}\big{]}=\operatorname{Tr}\big{[}(\sum_{k\in[K]}\lambda_{k}D_{k}^{-1})^{-1}\big{]}$
$\displaystyle=\sum_{j\in[n]}\frac{1}{\sum_{k\in[K]}d_{kj}^{-1}\lambda_{k}}.$

By exploiting Proposition 5.5, we have the reformulation of the feasible set

$\\{\lambda\in\Delta_{K}:\ell_{\mathcal{Z}}(\widehat{\mu}-\mu_{1})\leq\delta\\}=\left\\{\lambda\in\Delta_{K}:\displaystyle-\delta\leq\frac{\sum_{k\in[K]}c_{kj}\lambda_{k}}{\sum_{k\in[K]}d_{kj}^{-1}\lambda_{k}}-v_{j}^{\top}\mu_{1}\leq\delta\quad\forall
j\in[n]\right\\}.$

Combine the above formulae, problem (A.1) becomes

$\begin{array}[]{cl}\min&\displaystyle\sum_{j\in[n]}\frac{1}{\sum_{k\in[K]}d_{kj}^{-1}\lambda_{k}}\\\
\operatorname{s.t.}&\lambda\in\Delta_{K},\\\
&\displaystyle-\delta\leq\frac{\sum_{k\in[K]}c_{kj}\lambda_{k}}{\sum_{k\in[K]}d_{kj}^{-1}\lambda_{k}}-v_{j}^{\top}\mu_{1}\leq\delta\quad\forall
j\in[n],\end{array}$ (A.2)

Because each element $d_{kj}$ is strictly positive and any feasible solution
$\lambda$ should satisfy $\lambda\neq 0$, the fractional constraint of problem
(A.2) can be rewritten as systems of linear inequality constraints

$\forall
j\in[n]:\qquad\left\\{\begin{array}[]{l}\sum_{k\in[K]}c_{kj}\lambda_{k}\leq(\delta+v_{j}^{\top}\mu_{1})\sum_{k\in[K]}d_{kj}^{-1}\lambda_{k}\\\
\sum_{k\in[K]}c_{kj}\lambda_{k}\geq(v_{j}^{\top}\mu_{1}-\delta)\sum_{k\in[K]}d_{kj}^{-1}\lambda_{k}.\end{array}\right.$

Moreover, each fractional term in the objective function of (A.2) can be
reformulated as a second-order cone following Section 2.3 in Lobo,
Vandenberghe, Boyd, and Lebret (1998). This completes the proof. ∎

###### Proof of Lemma 5.7.

By exploiting the definition of $\mathcal{N}_{V}$, we can parametrize the
optimal solution as
$\pi_{k}^{\prime}\sim\mathcal{N}(\mu_{k}^{\prime},\Sigma_{k}^{\prime})$ where
$\Sigma_{k}^{\prime}=VD_{k}V^{\top}$ and $D_{k}=\operatorname{diag}(d_{k})$ is
a diagonal matrix.

First, consider Assertion (i). By utilizing the expression of the KL
divergence between Gaussian distributions in Definition 5.1, we have

$\begin{array}[]{rcl}(\mu_{k}^{\prime},\Sigma_{k}^{\prime},d_{k})=\arg&\min&(\widehat{\mu}_{k}-\mu_{k})^{\top}\widehat{\Sigma}_{k}^{-1}(\widehat{\mu}_{k}-\mu_{k})+\log\det\widehat{\Sigma}_{k}+\operatorname{Tr}\big{[}\widehat{\Sigma}_{k}^{-1}\Sigma_{k}\big{]}\\\
&\operatorname{s.t.}&\widehat{\mu}_{k}\in\mathbb{R}^{n},~{}\widehat{\Sigma}_{k}\in\mathbb{S}_{++}^{n},~{}\widehat{d}_{k}\in\mathbb{R}_{++}^{n},~{}\widehat{\Sigma}_{k}=V\operatorname{diag}(\widehat{d}_{k})V^{\top}.\end{array}$

For any $\widehat{\Sigma}_{k}\in\mathbb{S}_{++}^{n}$, the optimal solution in
the variable $\widehat{\mu}_{k}$ is $\mu_{k}^{\prime}=\mu_{k}$. The optimal
solution in the variable $\widehat{d}_{k}$ satisfies

$d_{k}=\arg\min\limits_{\widehat{d}_{k}\in\mathbb{R}_{++}^{n}}~{}\log\det\operatorname{diag}(\widehat{d}_{k})+\operatorname{Tr}\big{[}\operatorname{diag}(\widehat{d}_{k}^{-1})V^{\top}\Sigma_{k}V\big{]},$

which admits the optimal solution

$(d_{k})_{j}=v_{j}^{\top}\Sigma_{k}v_{j}\quad\forall j\in[n].$

Next, consider Assertion (ii). By utilizing the expression of the KL
divergence between Gaussian distributions in Definition 5.1, we have

$\begin{array}[]{rcl}(\mu_{k}^{\prime},\Sigma_{k}^{\prime},d_{k})=\arg&\min&(\widehat{\mu}_{k}-\mu_{k})^{\top}(\Sigma_{k})^{-1}(\widehat{\mu}_{k}-\mu_{k})-\log\det\widehat{\Sigma}_{k}+\operatorname{Tr}\big{[}\widehat{\Sigma}_{k}(\Sigma_{k})^{-1}\big{]}\\\
&\operatorname{s.t.}&\widehat{\mu}_{k}\in\mathbb{R}^{n},~{}\widehat{\Sigma}_{k}\in\mathbb{S}_{++}^{n},~{}\widehat{d}_{k}\in\mathbb{R}_{++}^{n},~{}\widehat{\Sigma}_{k}=V\operatorname{diag}(\widehat{d}_{k})V^{\top}.\end{array}$

Because $\Sigma_{k}\in\mathbb{S}_{++}^{p}$, the optimal solution in the
variable $\widehat{\mu}_{k}$ is $\mu_{k}^{\prime}=\mu_{k}$. The optimal
solution in the variable $\widehat{d}_{k}$ satisfies

$d_{k}=\arg\min\limits_{\widehat{d}_{k}\in\mathbb{R}_{++}^{n}}~{}-\log\det\operatorname{diag}(\widehat{d}_{k})+\operatorname{Tr}\big{[}\operatorname{diag}(\widehat{d}_{k})V^{\top}(\Sigma_{k})^{-1}V\big{]},$

which admits the optimal solution

$(d_{k})_{j}=[V^{\top}(\Sigma_{k})^{-1}V]_{jj}\quad\forall j\in[n].$

This completes the proof. ∎

###### Proof of Proposition 6.5.

Note that by Proposition 6.4, the first moment of the solution $\pi^{\star}$
to (14) satisfy

$\mathbb{E}_{\pi^{\star}}[\theta]=\widehat{\mu}=\sum_{k\in[K]}\lambda_{k}\mu_{k}.$

Now $(i)$ and $(iii)$ follow from definition. For $(ii)$, note that

$\Lambda_{\ell_{\mathcal{Z}},\delta}^{W}(\\{\pi_{k}\\})=\left\\{\lambda\in\Delta_{K}:\sup_{Az\leq
b}z^{\top}\left(\sum_{k\in[K]}\lambda_{k}\mu_{k}-\mu_{1}\right)\leq\delta\right\\}.$

The dual linear programming of $\sup_{Az\leq
b}z^{\top}\left(\sum_{k\in[K]}\lambda_{k}\mu_{k}-\mu_{1}\right)$ is

$\begin{array}[]{cl}\inf&b^{\top}w\\\
\operatorname{s.t.}&A^{\top}w=\sum_{k\in[K]}\lambda_{k}\mu_{k}-\mu_{1}\\\
&w\geq 0,\end{array}$

and by strong duality, the primal and dual optimal objective values are equal.
Therefore we can reformulate $\sup_{Az\leq
b}z^{\top}\left(\sum_{k\in[K]}\lambda_{k}\mu_{k}-\mu_{1}\right)\leq\delta$ as
the existence of $w\geq 0$ such that $b^{\top}w\leq\delta$ and
$A^{\top}w=\sum_{k\in[K]}\lambda_{k}\mu_{k}-\mu_{1}$. ∎

###### Proof of Theorem 6.6.

Because the covariance matrices $\Sigma_{k}$ are simultaneously
diagonalizable, Remark 4.4 in Álvarez-Esteban, del Barrio, Cuesta-Albertos,
and Matrán (2016) implies that the covariance matrix $\widehat{\Sigma}$ is
known in closed form

$\widehat{\Sigma}=\big{(}\sum_{k\in[K]}\lambda_{k}\Sigma_{k}^{\frac{1}{2}}\big{)}^{2}.$

Since the trace operator is a linear operator, we find

$\operatorname{Tr}\big{[}\widehat{\Sigma}\big{]}=\sum_{i,j\in[K]}\lambda_{i}\lambda_{j}\operatorname{Tr}\big{[}\Sigma_{i}^{\frac{1}{2}}\Sigma_{j}^{\frac{1}{2}}\big{]}.$

This completes the proof. ∎

###### Proof Theorem 7.3.

By exploiting Proposition 7.2, the first and second moments of the posterior
distribution $\pi^{\star}=\mathcal{C}_{\lambda}^{B}(\\{\pi_{k}\\})$ satisfy

$\mathbb{E}_{\pi^{\star}}[\theta]=\sum_{k\in[K]}\lambda_{k}\mathbb{E}_{\pi_{k}}[\theta]=\sum_{k\in[K]}\lambda_{k}\mu_{k},$
(A.3a) and
$\mathbb{E}_{\pi^{\star}}[\theta\theta^{\top}]=\sum_{k\in[K]}\lambda_{k}\mathbb{E}_{\pi_{k}}[\theta\theta^{\top}]=\sum_{k\in[K]}\lambda_{k}(\mu_{k}\mu_{k}^{\top}+\Sigma_{k}).$
Thus, the covariance matrix under $\pi^{\star}$ is
$\mathrm{Cov}_{\pi^{\star}}[\theta]=\mathbb{E}_{\pi^{\star}}[\theta\theta^{\top}]-\mathbb{E}_{\pi^{\star}}[\theta]\mathbb{E}_{\pi^{\star}}[\theta]^{\top}=\sum_{k\in[K]}\lambda_{k}(\mu_{k}\mu_{k}^{\top}+\Sigma_{k})-(\sum_{k\in[K]}\lambda_{k}\mu_{k})(\sum_{k\in[K]}\lambda_{k}\mu_{k})^{\top}.$
(A.3b)

Recall that $\lambda^{\star}$ solves

$\begin{array}[]{cl}\min&\mathrm{Tr}(\mathrm{Cov}_{\pi^{\star}}[\theta])\\\
\operatorname{s.t.}&\lambda\in\Delta_{K}\\\
&\ell_{\mathcal{Z}}(\mathds{E}_{\pi^{\star}}[\theta]-\mu_{1})\leq\delta.\end{array}$

The reformulation (19) follows by substituting (A.3a) and (A.3b) into the
above optimization problem. ∎

##  Appendix B Simulation Study

Figure A.1: Look-ahead-bias and variance with full Wasserstein mechanism on
simulated data

(a) Missing completely at random

(b) Missing at random

(c) Block missing

(d) Missing by value

These figures show the look-ahead-bias measured by the regret ECBias2, the
variance ECVar and the expected conditional mean squared error ECMSE for
simulated data with the full Wasserstein mechanism for different weights
$\delta$, which control the trade-off between look-ahead-bias and variance.
Figure A.2: Look-ahead-bias and variance with restricted Wasserstein mechanism
on simulated data

(a) Missing completely at random

(b) Missing at random

(c) Block missing

(d) Missing by value

These figures show the look-ahead-bias measured by the regret ECBias2, the
variance ECVar and the expected conditional mean squared error ECMSE for
simulated data with the restricted Wasserstein mechanism for different weights
$\delta$, which control the trade-off between look-ahead-bias and variance.

##  Appendix C Empirical Results

###  C.1 Imputation of Individual Stock Returns

Figure A.3: Look-ahead-bias and variance with full Wasserstein mechanism for
individual stocks

(a) Missing completely at random

(b) Missing at random

(c) Block missing

(d) Missing by value

These figures show the look-ahead-bias measured by the regret ECBias2, the
variance ECVar and the expected conditional mean squared error ECMSE for
individual stocks with the full Wasserstein mechanism for different weights
$\delta$, which control the trade-off between look-ahead-bias and variance.
The sample are daily returns of 10 representative large capitalization stocks
from 01/01/2015 to 01/01/2017. Figure A.4: Look-ahead-bias and variance with
restricted Wasserstein mechanism for individual stocks

(a) Missing completely at random

(b) Missing at random

(c) Block missing

(d) Missing by value

These figures show the look-ahead-bias measured by the regret ECBias2, the
variance ECVar and the expected conditional mean squared error ECMSE for
individual stocks with the restricted Wasserstein mechanism for different
weights $\delta$, which control the trade-off between look-ahead-bias and
variance. The sample are daily returns of 10 representative large
capitalization stocks from 01/01/2015 to 01/01/2017.

###  C.2 Imputation of Industry Portfolio Returns

Figure A.5: Look-ahead-bias and variance with full Wasserstein mechanism for
industry portfolios

(a) Missing completely at random

(b) Missing at random

(c) Block missing

(d) Missing by value

These figures show the look-ahead-bias measured by the regret ECBias2, the
variance ECVar and the expected conditional mean squared error ECMSE for
industry portfolios with the full Wasserstein mechanism for different weights
$\delta$, which control the trade-off between look-ahead-bias and variance.
The sample are daily returns of 10 industry portfolios from 01/03/2017 to
01/02/2019. Figure A.6: Look-ahead-bias and variance with restricted
Wasserstein mechanism for industry portfolios

(a) Missing completely at random

(b) Missing at random

(c) Block missing

(d) Missing by value

These figures show the look-ahead-bias measured by the regret ECBias2, the
variance ECVar and the expected conditional mean squared error ECMSE for
industry portfolios with the restricted Wasserstein mechanism for different
weights $\delta$, which control the trade-off between look-ahead-bias and
variance. The sample are daily returns of 10 industry portfolios from
01/03/2017 to 01/02/2019.