# Time-Series Imputation with Wasserstein Interpolation for Optimal Look-
Ahead-Bias and Variance Tradeoff

Jose Blanchet, Fernando Hernandez, Viet Anh Nguyen, Markus Pelger, Xuhui Zhang
Department of Management Science and Engineering, Stanford University

###### Abstract.

Missing time-series data is a prevalent practical problem. Imputation methods
in time-series data often are applied to the full panel data with the purpose
of training a model for a downstream out-of-sample task. For example, in
finance, imputation of missing returns may be applied prior to training a
portfolio optimization model. Unfortunately, this practice may result in a
look-ahead-bias in the future performance on 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 posterior consensus distribution which optimally controls
the variance and look-ahead-bias trade-off in the imputation. We demonstrate
the benefit of our methodology both in synthetic and real financial data.

##  1\. Introduction

Missing data is a common practical problem often encountered in data
collection and, as such, is prevalent in virtually any area touched by data
science. There is a large body of literature which deals with general
imputation tools, see for example, [23]. Understanding the implications of
imputing time-series data in downstream out-of-sample tasks, however, is
significantly less studied. This is precisely the goal of this paper.

We propose a novel Bayesian approach which interpolates between the use of
globally observed data (i.e. data collected during the whole 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 partial data to impute).

Naturally, the philosophy applied in virtually every imputation method
involves extracting as much information as possible from the globally observed
data to infer the missing data. However, in time-series settings (and others
in which a downstream task occurs relative to some ordering), this philosophy
poses a problem. This is of particular importance in financial applications
with portfolio optimization as a leading example. Portfolio selection
procedures are designed to exploit subtle signals that can lead to favorable
future investment outcomes. As spurious signal may induce a poor strategy, the
investments are typically evaluated out-of-sample to avoid potential
overfitting by “looking into the future”.

However, missing data is a prevalent problem in financial data sets: A large
number of instruments is traded over the counter or infrequently and many
assets might only temporarily trade on public exchanges. If missing data is
imputed globally using all information, an out-of-sample portfolio selection
procedure on the imputed data will inherently suffer from a look-ahead-bias.
This look-ahead-bias can be avoided by only using past data for imputation and
hence avoid “looking into the future”. This means that sequentially more data
is added to imputation procedure. While this approach avoids the look-ahead-
bias problem, it ignores valuable information that could be used to improve
the error in the imputation procedure. While overfitting by using future
information has been studied in the finance literature [12, 3, 15], none of
the existing results so far explores how to formulate and quantify look-ahead-
bias and optimally trade off between look-ahead-bias and variance in
imputations for future out-of-sample tasks.

Our conceptual idea extends to general imputation in time-series that are
followed by a downstream task. We choose the finance setting as a particularly
relevant application. While a look-ahead-bias in the financial context is
generally understood as using data that is available only in the future, as
for example in [3], the exact quantification and definition is a challenging
problem. We adopt a Bayesian imputation framework in order to precisely define
look-ahead-bias as deviation from a specific baseline posterior distribution.
We assume a Bayesian model with the observed entries as realizations. Such a
model is endowed with a prior distribution for the underlying parameters. The
full Bayesian model generates a corresponding posterior both for the
parameters and the missing observations conditional on the observed data. This
type of formulation is standard in the multiple imputation literature [28,
23]. The multiple imputation framework, from which we borrow inspiration,
imputes by sampling (multiple times) from the posterior, creating multiple
versions of imputed data set, each of which can be used for relevant
downstream tasks.

The 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 the Bayesian approach. Each of them corresponds to a
different posterior. The posterior computed by fitting the Bayesian model only
on the training set (which we call baseline) is assumed to be look-ahead-bias-
free, and the look-ahead-bias incurred by any other posterior is defined as
the bias relative to the baseline posterior. Our approach differs from the
standard multiple imputation literature in that we are interested in combining
these two natural posteriors to minimize the bias and variance trade-off in
the estimation of the returns. We therefore will obtain a “consensus”
posterior that is optimal in the sense of this trade-off. Most of the
literature in multiple imputation, in fact, would naturally advocate using the
full panel data.

We illustrate our novel conceptual framework in an i.i.d. Gaussian model with
a known covariance matrix. Thus, the moments of the posterior of missing
entries reduce to the moments of the posterior of the mean of the Gaussian
model. We use a variational approach to synthesize the two posteriors (fitted
from the training set and the full panel) of the mean parameter. A single
consensus posterior is selected by computing a Wasserstein barycenter
interpolation of the two posteriors. The barycenter interpolation introduces
weights which prescribe the importance attached to the mutual information
between the consensus and each of the posteriors. This set of weights is
further optimized to find the optimal consensus posterior according to a
criterion that minimizes the variance of the consensus posterior while
constraining its look-ahead-bias not to exceed some threshold. The threshold
parameter can then be either chosen by cross validation or by minimizing the
estimated mean squared error over all consensus posteriors as a function of
the threshold parameter.

Finally, we empirically validate the benefit of our methodology with the
downstream task of portfolio optimization, in which 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 natural posteriors. We demonstrate
our method on both synthetic and real financial data, and under different
practically relevant missing patterns, which include the simple case of
missing at random, a dependent block missing structure, and when the missing
pattern depends on past realizations of returns.

Contributions. In summary, our contributions are

  * •

We propose a novel conceptual framework to formalize and quantify the look-
ahead-bias in the context of data imputation. We combine the different layers
of future information used in the imputation via a Bayesian consensus (between
posterior distributions) in an optimal way that trades off between the look-
ahead-bias and variance.

  * •

We derive a tractable two-step optimization procedure to find the optimal
consensus posterior using the Wasserstein barycenter interpolation. The first
step finds the best possible distribution to summarize the different posterior
models for a given set of weights. The second step optimizes these weights
according to the bias-variance criterion that minimizes the variance while
explicitly controlling the look-ahead-bias.

##  2\. Related Literature

We briefly review the current imputation methods. Broadly speaking, the
various imputation methods can be classified into “single imputation” versus
“multiple imputation”. Single imputation refers to the use of a specific
number (i.e., a best guess) in place of each missing value, which includes for
example, k-nearest neighbor, mean/median imputation, smoothing, interpolation
[21] and spline. Matrix completion and matrix factorization [5, 6, 7, 24, 20]
are often used to impute panel data, as well as latent factor modelling [10,
2, 33]. Recently, there is a surge in the application of the recurrent neural
network for imputation [9, 11, 22, 8], as well as the generative adversarial
network [34].

Multiple imputation generalizes the single imputation procedures in that the
missing entries are filled-in with multiple guesses instead of one single best
guess or estimate, accounting for the uncertainty involved in imputation.
Multiple imputation is first developed for non-response in survey sampling
[28]. Since then it has been extended to time-series data [17, 16], where a
key new element is to preserve longitudinal consistency in imputation. [17]
uses smooth basis functions to increase smoothness of imputation. [16] uses
gap-filling sequence of imputation for categorical time-series data and
produces smooth patterns of transitions. The flexibility and predictive
performance of multiple imputation method have been successfully demonstrated
in a variety of data sets, such as Industry and Occupation Codes [13], GDP in
African nation [17], and concentrations of pollutants in the arctic [18].

The theory of multiple imputation is based on the ignorability assumption,
i.e., the data is missing at random. [14] found that earnings data is not
missing at random, and hence there is a bias (which they call attrition bias)
in estimates of earnings properties. [4] proposed an approach when the random
missing pattern depends on observable variables. [19, 26, 31] use logistic
regression models to model the missing data distribution.

Notation. 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$.

##  3\. The Bias-Variance Multiple Imputation Framework

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, for all $i\in[n]:$

$(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 projects111 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.$
(1)

where $\theta,\varepsilon_{t}\in\mathbb{R}^{n}$ and $M_{t}\in\\{0,1\\}^{n}$.
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. The panel data $\mathcal{D}$ with missing entries (red). 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 $\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 versions of the data set. We make the following assumptions.

######  Assumption 3.1.

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

######  Assumption 3.2.

The missing pattern satisfies the ignorability assumption [23], 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}})~{}\forall
t\in[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$.222Using the construction of $P_{M_{t}}$ in
Footnote 1, 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
$\mathcal{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. 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. For convenience we
denote $T_{1}=T^{train}$ and $T_{2}=T$. The data set truncated to time
$T_{k},~{}k=1,2$ is denoted as
$\mathcal{D}_{k}=\\{(X_{t},M_{t}),t\in[T_{k}]\\}$, so that $\mathcal{D}_{1}$
coincides with the training data set while $\mathcal{D}_{2}$ coincides with
the whole data set $\mathcal{D}$. The $k$-th posterior for $\theta$, denoted
by $\pi_{k}$, is formulated conditional on $\mathcal{D}_{k}$. The following
theorem shows that $\pi_{k}$ are Gaussians under Assumptions 3.1 to 3.3. All
proofs in this paper are relegated to the supplemental material.

######  Theorem 3.4.

Under Assumptions 3.1 to 3.3, the posterior distribution of
$\theta|\mathcal{D}_{k}$ is governed by
$\pi_{k}\sim\mathcal{N}_{n}(\mu_{k},\Sigma_{k}),~{}k=1,2$, where if $\pi_{0}$
is non-informative, then

$\displaystyle\Sigma_{k}$
$\displaystyle=\left(\sum_{t\in[T_{k}]}(\mathcal{P}_{M_{t}}^{\perp})^{-1}(\Omega_{X_{t}}^{-1})\right)^{-1},$
$\displaystyle\mu_{k}$
$\displaystyle=\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),$ (2)

and if $\pi_{0}$ is $\mathcal{N}_{n}(\mu_{0},\Sigma_{0})$, then

$\displaystyle\Sigma_{k}\\!$
$\displaystyle=\left(\Sigma_{0}^{-1}+\sum_{t\in[T_{k}]}(\mathcal{P}_{M_{t}}^{\perp})^{-1}(\Omega_{X_{t}}^{-1})\right)^{-1},$
$\displaystyle\mu_{k}\\!$
$\displaystyle=\\!\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 zeros333This
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})$..

The product of the posterior generator module $(G)$ is two elementary
posteriors $\\{\pi_{k}\\}_{k\in[2]}$, where $\pi_{k}$ is the posterior of
$\theta|\mathcal{D}_{k}$. The definition of the set $\mathcal{D}_{1}$ and
$\mathcal{D}_{2}$ implies that the posterior $\pi_{2}$ carries look-ahead-bias
relative to the posterior $\pi_{1}$, meanwhile also has less variance. The two
posteriors are transmitted to the consensus mechanism $(C)$ to form an
aggregated posterior $\pi^{\star}$ that strikes a balance between the look-
ahead-bias and variance. 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$.

######  Proposition 3.5.

Under Assumption 3.3, the distribution of $Y_{t}$ conditional on $\theta$ and
$(X_{t},M_{t})_{t\in[T]}$ is independent across $t$, and

$Y_{t}|\theta,(X_{t},M_{t})_{t\in[T]}\sim\mathcal{N}_{\mathrm{dim}(Y_{t})}\left(\tilde{\theta}_{t},\tilde{\Omega}_{t}\right),$
(4)

where
$\tilde{\theta}_{t}=\theta_{Y_{t}}+\Omega_{Y_{t}X_{t}}\Omega_{X_{t}}^{-1}(X_{t}-\theta_{X_{t}})$
and
$\tilde{\Omega}_{t}=\Omega_{Y_{t}}-\Omega_{Y_{t}X_{t}}\Omega_{X_{t}}^{-1}\Omega_{X_{t}Y_{t}}$.

Proposition 3.5 indicates 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},\theta\sim\pi^{\star},\eta_{t}\sim\mathcal{N}_{\mathrm{dim}(Y_{t})}(0,S_{t}),\theta\perp\\!\\!\\!\perp\eta_{t},$
(5)

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 (4).

The variability in the sampling of $\xi_{t}$ using (5) 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 in Theorem 3.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
look-ahead-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 1.

Algorithm 1 Conditional Expectation Bayesian Imputation

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

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

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

Compute $A_{t}$, $b_{t}$ for model (5) using $(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 two posteriors
$\pi_{1},\pi_{2}$ from the posterior generator module (G), then synthesizes a
unique posterior $\pi^{\star}$ and transmits $\pi^{\star}$ to the sampler
module (S). We provide a formal definition of a consensus mechanism. First, we
define the Wasserstein distance.

######  Definition 4.1 (Wasserstein distance).

For $p\geq 1$ the $p$-Wasserstein distance for two probability measures
$\pi_{1},\pi_{2}\in\mathcal{M}$ is

$W_{p}(\pi_{1},\pi_{2})=\left(\inf_{\pi\in\Pi(\pi_{1},\pi_{2})}\int_{\mathbb{R}^{2n}}\|x-y\|_{2}^{p}\mathrm{d}\pi(x,y)\right)^{1/p}.$

where $\Pi(\pi_{1},\pi_{2})$ is the set of all probability measures on
$\mathbb{R}^{2n}$ that have marginals $\pi_{1}$ and $\pi_{2}$.

We focus on the $2$-Wasserstein distance due to its analytical tractability.
Utilizing this distance, we define the Wasserstein consensus as follows.

######  Definition 4.2 (Consensus mechanism).

The Wasserstein consensus mechanism
$\mathcal{C}_{\lambda}:\mathcal{M}^{2}\to\mathcal{M}$ parametrized by
$\lambda$ outputs a unifying posterior $\pi^{\star}$ as

$\mathcal{C}_{\lambda}(\\{\pi_{k}\\})=\arg\min\limits_{\pi\in\mathcal{M}}\lambda_{1}W_{2}^{2}(\pi,\pi_{1})+\lambda_{2}W_{2}^{2}(\pi,\pi_{2}).$
(6)

The space of all possible consensus mechanism induced by the simplex
$\Delta=\\{\lambda\in\mathbb{R}_{+}^{2}:\lambda_{1}+\lambda_{2}=1\\}$ is
represented by
$\mathbb{C}\triangleq\\{\mathcal{C}_{\lambda}:\lambda\in\Delta\\}$.
Wasserstein distance enjoyes many sound properties, both in theory [32] and in
practice [29, 27]. Problem (6) computes the probability measure that minimizes
the sum of the distances to a set of target measures. The optimal solution is
usually called the Wasserstein barycenter [1], or McCann’s interpolation [25]
in the case of two target measures. The parameter vector $\lambda$ can be
thought of as the vector of weights assigning to the set of posteriors
$\\{\pi_{k}\\}$. Problem (6) 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, it
is well documented that $\mathcal{C}_{\lambda}$ can be characterized
analytically [25, Example 1.7].

######  Proposition 4.3.

Suppose that $\pi_{k}\sim\mathcal{N}(\mu_{k},\Sigma_{k})$ with
$\Sigma_{k}\succ 0$ for $k=1,2$. For any $\lambda\in\Delta$, we have
$\mathcal{C}_{\lambda}(\\{\pi_{k}\\})=\mathcal{N}(\widehat{\mu},\widehat{\Sigma})$
with

$\displaystyle\widehat{\mu}$
$\displaystyle=\lambda_{1}\mu_{1}+\lambda_{2}\mu_{2}\in\mathbb{R}^{n},$
$\displaystyle\widehat{\Sigma}$
$\displaystyle=\left(\lambda_{1}I_{n}+\lambda_{2}\Phi\right)\Sigma_{1}\left(\lambda_{1}I_{n}+\lambda_{2}\Phi\right)\in\mathbb{S}_{++}^{n},$

where
$\Phi=\Sigma_{2}^{1/2}\left(\Sigma_{2}^{1/2}\Sigma_{1}\Sigma_{2}^{1/2}\right)^{-1/2}\Sigma_{2}^{1/2}$.

Next we identify the optimal consensus mechanism
$\mathcal{C}_{\lambda^{\star}}\in\mathbb{C}$ according to the notion of bias-
variance trade-off, which is ubiquitously used in (Bayesian) statistics. The
first step towards this goal is to quantify the variance of the aggregated
posterior. Let
$\mathds{V}\mathrm{ar}_{\mathcal{C}_{\lambda}(\\{\pi_{k}\\})}[\theta]\in\mathbb{R}^{n}$
be the vector of variances 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 individual
variances.

######  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

$\mathbbm{1}^{\top}\mathds{V}\mathrm{ar}_{\pi^{\star}}[\theta]\leq\mathbbm{1}^{\top}\mathds{V}\mathrm{ar}_{\mathcal{C}_{\lambda}(\\{\pi_{k}\\})}[\theta]\quad\forall\mathcal{C}_{\lambda}\in\mathbb{C},$

where $\mathbbm{1}$ is the vector of ones.

The second steps entails quantifying the bias of an aggregated posterior
$\pi^{\star}$. Ideally, this bias quantity is measured by

$\|\mathds{E}_{\pi^{\star}}[\theta]-\theta^{true}\|_{2},$ (7)

where $\theta^{true}$ is the true mean of the returns vector. Unfortunately,
the value $\theta^{true}$ is not available in practice, and hence quantifying
the true bias of an aggregated posterior is impossible. As an alternative, we
resort to a proxy which is obtained through the following more general
definition of the bias notion.

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

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

The tolerance parameter $\delta$ indicates how much of bias is accepted in the
aggregated posterior $\pi^{\star}$. 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 future information, 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. As a result we have the following notion of an
optimal consensus mechanism.

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

For a given input $\\{\pi_{k}\\}$. The consensus mechanism
$\mathcal{C}_{\lambda^{\star}}\in\mathbb{C}$ is optimal for $\\{\pi_{k}\\}$ if
$\lambda^{\star}$ solves

$\begin{array}[]{cl}\min&\mathbbm{1}^{\top}\mathds{V}\mathrm{ar}_{\mathcal{C}_{\lambda}(\\{\pi_{k}\\})}[\theta]\\\
\operatorname{s.t.}&\lambda\in\Delta\\\
&\|\mathds{E}_{\mathcal{C}_{\lambda}(\\{\pi_{k}\\})}[\theta]-\mathds{E}_{\pi_{1}}[\theta]\|_{2}\leq\delta.\end{array}$
(8)

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 look-ahead-bias.
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. It is important to
notice the sharp distinction between our consensus optimization problem (8)
and the existing literature. While the existing literature typically focuses
on studying the barycenter for a fixed weight parameter $\lambda$, problem (8)
searches for the best value of $\lambda$ that meets our optimality criteria.

Proposition 4.3 provides us with a tractable reformulation of problem (8),
which is the main result of this section.

######  Theorem 4.7.

Suppose that $\pi_{k}\sim\mathcal{N}(\mu_{k},\Sigma_{k})$ with
$\Sigma_{k}\succ 0$ for $k=1,2$. Let $\lambda^{\star}$ be the optimal solution
in the variable $\lambda$ of the following quadratic programming

$\begin{array}[]{cl}\min&\mathrm{Tr}(\Sigma_{1})\lambda_{1}^{2}+\mathrm{Tr}(\Sigma_{2})\lambda_{2}^{2}+2\mathrm{Tr}(\Sigma_{1}\Phi)\lambda_{1}\lambda_{2}\\\
\operatorname{s.t.}&\lambda\in\Delta\\\
&\lambda_{2}\|\mu_{1}-\mu_{2}\|_{2}\leq\delta,\end{array}$ (9)

where
$\Phi=\Sigma_{2}^{1/2}\left(\Sigma_{2}^{1/2}\Sigma_{1}\Sigma_{2}^{1/2}\right)^{-1/2}\Sigma_{2}^{1/2}$.
Then $\mathcal{C}_{\lambda^{\star}}$ is the optimal consensus mechanism in
$\mathbb{C}$.

##  5\. Numerical Experiment

In this section we design and carry out numerical experiments to validate the
look-ahead-bias and variance trade-off that is captured by our proposed BVMI
framework. All optimization problems are modeled in Python and solved by the
solver GUROBI on an Intel Core i5 CPU (1.40 GHz) computer.

The panel of asset returns data (either synthetic or real) is separated into
three portions: training, testing, and an additional out-of-sample testing
periods, each having a length $T^{train},T^{test}$ and $T^{oos-test}$
respectively. In both synthetic and real data experiments we will manually
create missingess by applying masks to entries in the panel following certain
missing mechanisms. We only apply the masking procedure to the training
period, while we preserve the completeness of testing and out-of-sample
testing periods for evaluation purposes.

Following the BVMI framework, we feed the data set with missing values to the
posterior generator module, as we have detailed in Section 3.1. Recall that we
have a total of two posterior distributions, which are fitted using the
training period and training plus testing periods respectively. 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
synthetic data experiment, or fitted from the data in the case of real data
experiment. The generated multiple posterior distributions are then fed to the
consensus mechanism module to output an optimal consensus posterior, the
details of which we have provided in Sections 4. As we have remarked before,
the tolerance parameter $\delta$ in (8) (or equivalently in (9)) indicates how
much of look-ahead-bias is accepted in the aggregated posterior $\pi^{\star}$.
We experiment with the $\delta$ parameter selected to be $10$ values equi-
spaced between $0$ and $\delta_{max}$, where $\delta=0$ gives us the
aggregated posterior $\pi^{\star}$ matching the posterior $\pi_{1}$, and
$\delta=\delta_{max}$ gives us the aggregated posterior $\pi^{\star}$ matching
the posterior $\pi_{2}$. We can compute $\delta_{max}$ by
$\|\mu_{1}-\mu_{2}\|_{2}$. Corresponding to each selected $\delta$ parameter,
the aggregated posterior $\pi^{\star}$ is then output to the multiple
imputation sampler module, where we follow Algorithm 1 to create multiply-
imputed data sets.

The downstream task to be performed on the imputed data set is that of
portfolio optimization. We consider a greedy strategy to construct the
portfolio weights, that is, we calculate the $l_{2}$-normalized average
returns of the training period, and apply the weights on the testing and out-
of-sample testing periods. The portfolio weights $w\in\mathbb{R}^{n}$ is
chosen as

$w=\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}}=\frac{\sum_{t\in[T^{train}]}\tilde{Z}_{t}}{\|\sum_{t\in[T^{train}]}\tilde{Z}_{t}\|_{2}},$

where $\tilde{Z}_{t}$ denotes the imputed value of $Z_{t}$. When the portfolio
weights $w$ is applied on the testing period, we have the averaged daily
portfolio return

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

Similarly, when the portfolio weights $w$ is applied on the out-of-sample
testing period, we have the averaged daily portfolio return

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

The difference in averaged daily portfolio return between the two periods, or
regret, is defined as

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

Conditional on $\mathcal{D}$ (the realization of observed data and missing
masks), $\Delta R$ is a random variable where the randomness comes from
multiple imputations. The randomness in the conditional moments
$\mathbb{E}[\Delta R|\mathcal{D}]$ and $\mathbb{V}\mathrm{ar}[\Delta
R|\mathcal{D}]$ 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}]].$

Note that we consider one-sided bias in the above formulation. We call the
first term above as ECBias2, and the second term above as ECVar.

In the following, we make further specifications of the experiments depending
on whether we work with synthetic or real data. We present the (qualitative)
trend of ECMSE, ECBias2 and ECVar as $\delta$ varies, under three different
missing mechanisms which include missing at random, a dependent block missing
structure, and when the missing pattern depends on realizations of returns.
The results show that our proposed BVMI framework has the potential to
outperform naive posteriors in terms of ECMSE.

(a) Missing completely at random

(b) Block missing

(c) Missing by value

Figure 3. Experimental results on synthetic data set

(a) Missing completely at random

(b) Block missing

(c) Missing by value

Figure 4. Experimental results on real data set

###  5.1. Experiments with Synthetic Data

We consider number of stocks $n=100$, and horizons
$T^{train}=T^{test}=100,T^{oos-test}=1000$. We generate the returns vector
from a factor model, namely, the vectors
$Z_{t},t=1,\ldots,T^{train}+T^{test}+T^{oos-test}$ are generated i.i.d. from
$\mathcal{N}_{n}(\theta,\Omega)$, where for $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 three 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 $70\%$; block missing \- we mask the first $60\%$ 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.25$.

For each missing mechanism, we repeat $100$ independent experiments, where in
each experiment the sampled panel is multiply-imputed $100$ times, and we
carry out the downstream task on each imputed data set. The trends of ECMSE,
ECBias2 and ECVar with respect to the $\delta$ parameter are summarized in
Figure 3. We see that in all three cases of missing mechanisms, the bias
component ECBias2 is increasing as $\delta/\delta_{max}$ transitions from $0$
to $1$, while the variance component ECVar is decreasing. For the first two
cases, the sum of the two components ECMSE achieves a minimum for some non-
trivial selection of $\delta$, showing that our method outperforms the naive
posteriors (i.e., the posteriors corresponding to $\delta=0$ or
$\delta=\delta_{\max}$).

###  5.2. Experiments with Real Data

We choose the data as the $n=49$ industry portfolio returns downloaded from
Kenneth French’s website 444Retrieved on Jan. 24th. 2021 from
https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html in a
duration from Jan. 3nd. 2017 to Dec. 31st. 2018. We consider $T^{train}=200$,
$T^{test}=100$ and $T^{oos-test}=100$. The covariance matrix $\Omega$ in
Assumption 3.3 is fitted using the complete panel (a total of $502$ days). As
before, we consider three types of missing mechanisms. However, we use the
data set differently in each case.

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

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

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

The trends of ECMSE, ECBias2 and ECVar with respect to the $\delta$ parameter
are summarized in Figure 4. We see that in all three cases of missing
mechanisms, the bias component ECBias2 is increasing as $\delta/\delta_{max}$
transitions from $0$ to $1$, while the variance component ECVar is decreasing.
The sum of the two components ECMSE achieves a minimum for some non-trivial
selection of $\delta$ for the last two cases, showing that our method has the
potential to outperform the naive posteriors (i.e., the posteriors
corresponding to $\delta=0$ or $\delta=\delta_{max}$).

#### Acknowledgments.

We gratefully acknowledge the financial support and feedback of MSCI and its
data science team.

## References

  * [1] M. Agueh and G. Carlier, Barycenters in the wasserstein space, SIAM Journal on Mathematical Analysis, 43 (2011), pp. 904–924.
  * [2] J. Bai and S. Ng, Matrix completion, counterfactuals, and factor analysis of missing data, arXiv:1910.06677, (2020).
  * [3] R. W. Banz and W. J. Breen, Sample-dependent results using accounting and market data: Some evidence, The Journal of Finance, 41 (1986), pp. 779–793.
  * [4] C. Breunig, M. Kummer, J. Ohnemus, and S. Viete, Information technology outsourcing and firm productivity: eliminating bias from selective missingness in the dependent variable, The Econometrics Journal, 23 (2019), pp. 88–114.
  * [5] J.-F. Cai, E. J. Candès, and Z. Shen, A singular value thresholding algorithm for matrix completion, SIAM Journal on Optimization, 20 (2010), pp. 1956–1982.
  * [6] E. J. Candès and B. Recht, Exact matrix completion via convex optimization, Foundations of Computational Mathematics, 9 (2009), p. 717.
  * [7] E. J. Candes and T. Tao, The power of convex relaxation: Near-optimal matrix completion, IEEE Transactions on Information Theory, 56 (2010), pp. 2053–2080.
  * [8] W. Cao, D. Wang, J. Li, H. Zhou, Y. Li, and L. Li, Brits: Bidirectional recurrent imputation for time series, in Proceedings of the 32nd International Conference on Neural Information Processing Systems, 2018, p. 6776–6786.
  * [9] Z. Che, S. Purushotham, K. Cho, D. Sontag, and Y. Liu, Recurrent neural networks for multivariate time series with missing values, Nature Scientific Reports, 8 (2018), p. 6085.
  * [10] Y. Chen, J. Fan, C. Ma, and Y. Yan, Inference and uncertainty quantification for noisy matrix completion, Proceedings of the National Academy of Sciences, 116 (2019), pp. 22931–22937.
  * [11] E. Choi, M. T. Bahadori, A. Schuetz, W. F. Stewart, and J. Sun, Doctor ai: Predicting clinical events via recurrent neural networks, vol. 56 of Proceedings of Machine Learning Research, 2016, pp. 301–318.
  * [12] J. Choi, C. Liu, and J. Woo, An efficient approach for removing look-ahead bias in the least square monte carlo algorithm: Leave-one-out, arXiv:1810.02071, (2018).
  * [13] C. C. Clogg, D. B. Rubin, N. Schenker, B. Schultz, and L. Weidman, Multiple imputation of industry and occupation codes in census public-use samples using bayesian logistic regression, Journal of the American Statistical Association, 86 (1991), pp. 68–78.
  * [14] P. Easton, M. Kapons, P. Kelly, and A. Neuhierl, Attrition bias and inferences regarding earnings properties; evidence from compustat data, SSRN, (2020).
  * [15] H. Guo, H. Kassa, and M. F. Ferguson, On the relation between EGARCH idiosyncratic volatility and expected stock returns, Journal of Financial and Quantitative Analysis, 49 (2014), p. 271–296.
  * [16] B. Halpin, Multiple imputation for categorical time series, The Stata Journal, 16 (2016), pp. 590–612.
  * [17] J. Honaker and G. King, What to do about missing values in time series cross-section data, American Journal of Political Science, 54 (2010), pp. 561–581.
  * [18] P. K. Hopke, C. Liu, and D. B. Rubin, Multiple imputation for multivariate data with missing and below-threshold measurements: Time-series concentrations of pollutants in the arctic, Biometrics, 57 (2001), pp. 22–33.
  * [19] J. G. Ibrahim, S. R. Lipsitz, and M.-H. Chen, Missing covariates in generalized linear models when the missing data mechanism is non-ignorable, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61 (1999), pp. 173–190.
  * [20] R. H. Keshavan, A. Montanari, and S. Oh, Matrix completion from a few entries, IEEE Transactions on Information Theory, 56 (2010), pp. 2980–2998.
  * [21] D. M. Kreindler and C. J. Lumsden, The effects of the irregular sample and missing data in time series analysis, Nonlinear Dynamics Psychol Life Sci, 10 (2006), pp. 187–214.
  * [22] Z. C. Lipton, D. Kale, and R. Wetzel, Directly modeling missing data in sequences with rnns: Improved classification of clinical time series, vol. 56 of Proceedings of Machine Learning Research, 2016, pp. 253–270.
  * [23] R. Little and D. Rubin, Statistical analysis with missing data, Wiley, 2002.
  * [24] R. Mazumder, T. Hastie, and R. Tibshirani, Spectral regularization algorithms for learning large incomplete matrices, Journal of Machine Learning Research, 11 (2010), pp. 2287–2322.
  * [25] R. J. McCann, A convexity principle for interacting gases, Advances in Mathematics, 128 (1997), pp. 153–179.
  * [26] K. Morikawa, J. K. Kim, and Y. Kano, Semiparametric maximum likelihood estimation with data missing not at random, 45 (2017), pp. 393–409.
  * [27] O. Pele and M. Werman, Fast and robust earth mover’s distances, in International Conference on Computer Vision, vol. 9, 2009, pp. 460–467.
  * [28] D. B. Rubin, Multiple Imputation for Nonresponse in Surveys, Wiley, 1987\.
  * [29] Y. Rubner, L. Guibas, and C. Tomasi, The earth mover’s distance, multi-dimensional scaling, and color-based image retrieval, in Proceedings of the ARPA Image Understanding Workshop, 1997, pp. 661–668.
  * [30] M. J. Schervish, Theory of Statistics, Springer, 1995.
  * [31] A. Sportisse, C. Boyer, and J. Josse, Imputation and low-rank estimation with missing not at random data, Statistics and Computing, 30 (2020), pp. 1629–1643.
  * [32] C. Villani, Optimal Transport: Old and New, Springer Science & Business Media, 2008.
  * [33] R. Xiong and M. Pelger, Large dimensional latent factor modeling with missing observations and applications to causal inference, SSRN, (2020).
  * [34] J. Yoon, J. Jordon, and M. van der Schaar, GAIN: Missing data imputation using generative adversarial nets, vol. 80 of Proceedings of Machine Learning Research, 2018, pp. 5689–5698.

##  6\. Appendix

###### Proof of Theorem 3.4.

Let $f(\mathcal{D}_{k}|\theta)$ be the conditional density representing the
likelihood of observing $\mathcal{D}_{k}$ given the parameter $\theta$. Under
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}$. Now the Bayes theorem [30, Theorem 1.31] asserts that

$\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)}=\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)}$

where $\mathrm{d}\pi_{k}/\mathrm{d}\pi_{0}$ represents the Radon-Nikodym
derivative of $\pi_{k}$ with respect to $\pi_{0}$. Under Assumption 3.3
$X_{t}|\theta$ is Gaussian with mean $\theta_{X_{t}}$ and covariance matrix
$\Omega_{X_{t}}$. Thus the posterior $\pi_{k}$ satisfies

$\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}})$. Since
$\pi_{0}$ is either flat or Gaussian by Assumption 3.3, the posterior
$\pi_{k}$ also admits a density with respect to the Lebesgue measure on
$\mathbb{R}^{n}$. If $\pi_{0}$ follows $\mathcal{N}_{n}(\mu_{0},\Sigma_{0})$,
the density of $\pi_{k}$ 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, we see that
$\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),$

Since $\Omega_{X_{t}}^{-1}$ is positive definite, we have that
$(\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.

Similarly, if $\theta$ has a flat uninformative prior, then
$\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).$

We argue that the matrix inversion that defines $\Sigma_{k}$ is well-defined
under Assumption 3.1. 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 3.5.

Under Assumption 3.3, conditional on $\theta$, the joint distribution of
$X_{t},Y_{t}$ (i.e. $Z_{t}$) is Gaussian $\mathcal{N}_{n}(\theta,\Omega)$ and
is independent across $t$. The distribution of $Y_{t}$ conditional on $\theta$
and $(X_{t},M_{t})$ is thus a conditional Gaussian given by

$Y_{t}|\theta,X_{t},M_{t}=Y_{t}|\theta,(X_{t},M_{t})_{t\in[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).$

∎

###### Proof of Proposition 4.3.

From [1, Section 6.2] the Wasserstein barycenter for two target measures is
the well-known McCann’s interpolation [25]. If the two target measures are
Gaussians, then [25, Example 1.7] asserts the barycenter to be Gaussian with
the desired mean and covariance matrix. ∎

###### Proof of Theorem 4.7.

From Proposition 4.3, it is straightforward to compute

$\displaystyle\mathbbm{1}^{\top}\mathbb{V}ar_{\mathcal{C}_{\lambda}(\\{\pi_{k}\\})}[\theta]$
$\displaystyle=\mathrm{Tr}(\widehat{\Sigma})$
$\displaystyle=\mathrm{Tr}\left((\lambda_{1}I_{n}+\lambda_{2}\Phi)\Sigma_{1}(\lambda_{1}I_{n}+\lambda_{2}\Phi)\right)$
$\displaystyle=\lambda_{1}^{2}\mathrm{Tr}(\Sigma_{1})+\lambda_{1}\lambda_{2}(\mathrm{Tr}(\Sigma_{1}\Phi)+\mathrm{Tr}(\Phi\Sigma_{1}))+\lambda_{2}^{2}\mathrm{Tr}(\Phi\Sigma_{1}\Phi)$
$\displaystyle=\lambda_{1}^{2}\mathrm{Tr}(\Sigma_{1})+2\lambda_{1}\lambda_{2}\mathrm{Tr}(\Sigma_{1}\Phi)+\lambda_{2}^{2}\mathrm{Tr}(\Sigma_{2}),$

and

$\|\mathbb{E}_{\mathcal{C}_{\lambda}(\\{\pi_{k}\\})}[\theta]-\mathbb{E}_{\pi_{1}}[\theta]\|_{2}=\|\widehat{\mu}-\mu_{1}\|_{2}=\|(\lambda_{1}-1)\mu_{1}+\lambda_{2}\mu_{2}\|_{2}=\lambda_{2}\|\mu_{1}-\mu_{2}\|_{2}.$

Substitute the above expressions for problem (8) gives us the desired result.
∎