# Portfolio Optimization with Sparse Multivariate Modelling

PIER FRANCESCO PROCACCI 1 and TOMASO ASTE 1,2
1Department of Computer Science, UCL, Gower Street, London, WC1E6BT, UK
2Systemic Risk Centre, London School of Economics and Political Sciences,
London, WC2A 2AE, UK

(v1 released 28 Mar 2021)

###### Abstract

Portfolio optimization approaches inevitably rely on multivariate modeling of
markets and the economy. In this paper, we address three sources of error
related to the modeling of these complex systems: 1\. oversimplifying
hypothesis; 2\. uncertainties resulting from parameters’ sampling error; 3\.
intrinsic non-stationarity of these systems. For what concerns point 1. we
propose a $L_{0}$-norm sparse elliptical modeling and show that sparsification
is effective. The effects of points 2. and 3. are quantified by studying the
models’ likelihood in- and out-of-sample for parameters estimated over train
sets of different lengths. We show that models with larger off-sample
likelihoods lead to better performing portfolios up to when two to three years
of daily observations are included in the train set. For larger train sets, we
found that portfolio performances deteriorate and detaches from the models’
likelihood, highlighting the role of non-stationarity. We further investigate
this phenomenon by studying the out-of-sample likelihood of individual
observations showing that the system changes significantly through time.
Larger estimation windows lead to stable likelihood in the long run, but at
the cost of lower likelihood in the short-term: the ‘optimal’ fit in finance
needs to be defined in terms of the holding period. Lastly, we show that
sparse models outperform full-models in that they deliver higher out of sample
likelihood, lower realized portfolio volatility and improved portfolios’
stability, avoiding typical pitfalls of the Mean-Variance optimization.

###### keywords:

Portfolio Construction; Market States; Markowitz; Mean-Variance; Information
Filtering Networks; TMFG; Sparse inverse covariance; Correlation Structure.

C38, G61, G15, G17

##  1 Introduction

Quantitative approaches to asset management have accumulated unprecedented
popularity over the last few decades. Of all the algorithms and strategies
developed, portfolio selection models are among those that have received wider
attention. The essence of portfolio investing is to find the best way of
assigning weights to a given portfolio of assets to maximize future portfolio
returns while minimizing the investment risk. The exploration of this filed
starts with Markowitz’s Mean-Variance optimization process (Markowitz, 1952).
The theory of mean–variance-based portfolio selection is still today a
cornerstone of modern asset management. It rests on the presumption that
rational investors choose among risky assets purely on the basis of expected
return and risk, with risk measured as portfolio variance. The theoretical
foundation of this framework is sound if either: investors exhibit quadratic
utility, in which case they ignore non-normality in the data (Gollier, 2001),
or all the higher moments of the portfolio distribution can be expressed as a
function of mean and variance and hence all optimal solutions satisfy the
mean-variance criterion. Also, the“optimality” of the mean-variance portfolios
is based on the assumption that investors live in a one-period world, while in
reality they have an investment horizon that lasts longer than one period.
Markets indeed constantly change over time and investors are subject to
inflows/outflows forcing them to adjust their allocation and take corrective
actions.
Form a general, high level, perspective all portfolio optimization approaches
are based on a multivariate model of the variables in the market and the
economy. The optimization strategies are devised to maximize profits and
minimize risks based on such models. In modeling these complex systems there
are, however, several sources of inaccuracies and errors with the three main
ones being: 1\. oversimplifying hypothesis (such as the use of normal
distributions); 2\. uncertainties resulting from the estimation of the
parameters from datasets of limited sizes; 3\. intrinsic non-stationarity of
these systems, which makes in-sample estimations, based on past observations,
inadequate for the estimation of off-sample, future properties. Most likely
all three of these factors – and others – contribute to undermining the
predictive power of any attempt of modeling markets.

While a large deal of literature has been devoted to relaxing some of the most
unrealistic model assumptions (point 1.) the current main pitfall of portfolio
optimization is attributed to error maximization (point 2.). This effect has
long been established in literature (Michaud, 1989; Nawrocki, 1996).
Essentially, inputs into the mean-variance optimization are measured with
uncertainty, and the optimization procedure tends to pick those assets which
appear to have the most attractive features – but these are outlying cases
where estimation error is likely to be the highest, hence maximizing the
impact of estimation error on portfolios’ weights. The estimation error is
also amplified by market evolution which makes the training on the past not
fully representative of future market behavior (point 3.).

In this paper, we address all three sources of inaccuracies. For what concerns
point 1. we propose a $L_{0}$-norm topologically regularized sparse elliptical
modeling (Aste, 2020) and show that sparsification is effective. We quantify
the effects of estimation error and non-stationarity on portfolio performances
(point 2. and 3.) by assessing the goodness of models’ statistical likelihood
for estimates over train sets of different lengths. Specifically, we study how
the realized portfolio variance reacts to different out-of-sample likelihoods
of the input parameters and, particularly, to sparse models. Further, we
analyze how sparse precision matrices impact the magnitude and the stability
of portfolio weights.

The remainder of this paper is organized as follows: in Section 2 we briefly
review the theory around portfolio construction, highlighting the pitfalls on
assumptions and estimation error and the main solutions proposed in
literature; in Section 3 we outline our methodology and experiments design and
in Section 4 we present the results. Appendices 6 and 7 are devoted to
recalling some useful aspects of Elliptical distributions.

##  2 Literature Review

###  2.1 Modern Portfolio Theory

Considering a portfolio of $n$ assets with weights $\bm{W}=(w_{1},...,w_{n})$,
returns $\bm{R}=(\bm{r}_{1},...,\bm{r}_{n})$ and portfolio returns

$\bm{r}_{p}=\bm{WR}^{T},$ (1)

the standard Mean-Variance optimization problem consists in minimizing the
portfolios’ variance $\sigma_{p}$ for fixed levels of expected returns
$\mathbb{E}[\bm{r}_{p}]=\bar{r}_{p}$:

$\displaystyle\min_{\bm{W}}$ $\displaystyle\sigma_{p}^{2}=\bm{W\Sigma W}^{T}$
(2) s.t. $\displaystyle\mathbb{E}[\bm{r}_{p}]=\bar{r}_{p},$ and
$\displaystyle\bm{W}\mathbbm{1}=1,$

where $\bm{\Sigma}\in\mathbb{R}^{n\times n}$ is the assets’ covariance matrix
and $\mathbbm{1}\in\mathbb{R}^{n\times 1}$ is a basis column vector with all
elements equal 1. Solving for $\bm{W}$ for different values of $\bar{r}_{p}$,
one can obtain the optimal weights (i.e. the weights that minimize the
portfolio variance) corresponding to different portfolio expected returns
$\bar{r}_{p}$ yielding the so-called efficient frontier – i.e. the set of
optimal weights which provide the lowest variance for each level of expected
return.

As discussed in the introduction, this optimization is only concerned with the
first two moments of the distribution of portfolios’ returns and it does not
deal with multiperiod investment decisions. These pitfalls have largely been
discussed in literature. Kritzman (2000) provides a clear review of what are
the assumptions under which repeatedly investing in one-period-efficient
portfolios will also result in multiperiod efficiency. Also, many models have
been proposed to deal explicitly with multiperiod optimality (see, for
example, Li and Ng (2000) or Arditti and Levy (1976)). With respect to non-
normality of returns, a rich literature is available on both alternative
parametrization of data (Bamberg and Dorfleitner, 2001) and optimisation
frameworks that consider other distribution moments (Hogan and Warren, 1974;
Bawa, 1978) or other measures of risk/return (Harlow, 1991; Lwin et al.,
2017).

From the optimization problem in Eq.LABEL:min_variance it is also clear that
the optimization does not treat the error and uncertainty around the
parameters $\bm{\Sigma}$ and $\bm{\mu}$. The difference between the estimated
and true distribution parameters is called estimation error. It arises from
both the sampling procedure or availability of data and non-stationarity. The
error coming from sampling, also referred to as sampling error, is due to
parameters used in the portfolio optimization process being typically point
estimates – we can only expect these estimates to equal the true distribution
parameters if our sample is infinitely large. Assuming stationary data,
sampling error could be fixed by increasing the number of observations in the
estimation sample. Indeed, the convergence rate is in the inverse of the
square-root of the sample size, as dictated by the law of large numbers. This
would come handy in our times of increasing data availability. However, a
second source of estimation error comes from non-stationarity. A time series
is said to be non-stationary if its distribution parameters (or the
distribution itself) changes over time – in this case, extending the length of
observations might reduce the contribution of sampling error to estimation
error, but at the same time, it could increase that of non-stationarity
(Broadie, 1993).

Many techniques have been proposed in literature to deal with this phenomenon,
both relying on heuristic methods and decision-theoretic foundations (Scherer,
2007). Heuristic approaches mainly propose to constrain the optimization
problem in order to impose feasible optimal weights.

Michaud and Michaud (1998) address explicitly the sampling error proposing a
Monte Carlo based procedure called resampling. In order to model the
randomness of the input mean vector and covariance matrix, portfolio
resampling consists in repeatedly drawing from the return distribution given
by the point estimates and creating $n$ artificial new samples. For each
sample, an efficient frontier is estimated and the final, resampled-efficient
frontier is given by the average weight across all of the resampled
portfolios.

From a decision-theoretic perspective, Bayesian techniques have recently
played a primary role in literature. The rationale behind Bayesian statistics
for portfolio construction is to include non-sample information to tackle the
effect of parameter uncertainty on optimal portfolio choice. Instead of a
point estimate, Bayesian approaches produce a density function for the
parameters involved, by combining sample information (likelihood) with prior
belief, potentially coming from non-sample information. A special case of this
general approach is the seminal work of Black and Litterman (1992). In their
pioneering work, the authors assume assets’ returns to be normally distributed
with mean equal to the “equilibrium returns” (that is, the mean returns that
would return the market portfolio if used in a mean-variance optimization) and
combine this “sample” information with investors’ views on the assets. In this
way, in absence of an informative prior from investors, the model would return
the market or “equilibrium” portfolio. In presence of investors’ priors,
instead, the allocation would diverge from the equilibrium portfolio
accounting for investors’ views and proportionally to their confidence level.
Other than being highly appealing from a practitioner’s perspective, the model
proposed in Black and Litterman (1992) highlights the flexibility of the
Bayesian framework, with many sources of information that could potentially be
used in combination or to update the in-sample information. This is a very
active area of research with recant notable examples including (Scherer et
al., 2012) and (De Franco et al., 2019).

###  2.2 Precision matrix estimation with topological regularization

The literature concerning parameter estimation has greatly developed in recent
years. While the increasing availability of data has represented the fuel for
data-greedy machine learning models, it has also exacerbated problems related
to the curse of dimensionality and overfitting. In essence, the goal is to
extract the largest amount of information from data with a model that avoids
overfitting, generalizing well with new data. Furthermore, model
interpretability (or explainability) has become increasingly important and for
this purpose, simpler and sparser models with a smaller number of parameters
are preferred. To this extent, understanding the dependency structure among
variables has proven to be essential to capture the collective behaviour of
the systems and many techniques have been introduced in literature. Examples
are dimensionality reduction methods (Van Der Maaten et al., 2009) or
shrinkage techniques (Ledoit and Wolf, 2020; Friedman et al., 2008).

One possible approach to represent the set interactions in a complex system is
to model them as a network structure where the vertices are the system’s
elements and edges between vertices indicate the interactions between the
corresponding elements (Lauritzen, 1996). Information filtering networks
(Tumminello et al., 2005) aim at retrieving the relevant sub network of
interactions among the elements of the system. In the pioneering work of
Mantegna (1999), it was proposed to investigate financial systems by the
extraction of a minimal set of relevant interactions associated with the
strongest correlations belonging to the Minimum Spanning Tree (MST). The MST
structure is, however, a drastic filtering tool and is likely to discard
valuable information. Tumminello et al. (2005) and Aste and Di Matteo (2006)
expanded the concept by introducing the general idea of information filtering
networks. In particular, they show that graphs of different complexities can
be constructed by iteratively linking the most strongly connected nodes under
the constraint of generating planar or hyperbolic graphs. There is now a large
body of literature proving network filtering to be a powerful tool to
associate a sparse network to a high-dimensional dependency measure with
applications ranging from financial markets (Barfuss et al., 2016) to
biological systems (Song et al., 2012) and econophysics (Mantegna and Stanley,
2000). Recent developments extend the information filtering network
methodology to chordal graphs (Barfuss et al., 2016; Massara et al., 2015;
Massara and Aste, 2019). This allows to directly associate the information
filtering network with a positive definite sparse inverse covariance matrix
that is a $L_{0}$-norm topological regularization of the full covariance
estimate (Aste, 2020).

In our analysis, we use the TMFG-LoGo network filtering approach (Massara et
al., 2015; Barfuss et al., 2016). TMFG-LoGo approach has proven to be more
efficient and better performing, particularly when few data are available
(Barfuss et al., 2016; Aste and Di Matteo, 2017), with respect to $L_{1}$-norm
sparsification techniques such as GLASSO (Friedman et al., 2008). Moreover,
the TMFG procedure is computationally highly efficient and this is well suited
in our resampling experiment.

Despite the rich literature on portfolio construction techniques, to the best
of our knowledge, the link between a measure of estimation goodness and
portfolio performances is still not specifically treated in the literature.

##  3 Methodology

The estimation error is quantified by measuring how well the functional form
of the multivariate probability distribution, $f_{\bm{\theta}}(\bm{X})$
defined via the parameters $\bm{\theta}$ estimated in-sample, describes the
actual data out-of-sample. The statistical measure that describes how likely
are observations to belong to the estimated probability function is the
likelihood. The likelihood principle is a cornerstone of statistical analysis
in that maximum likelihood estimators are guaranteed to be asymptotically
efficient under mild conditions (Wald (1949), Cramér (1946), Daniels (1965)).
In this section, we introduce our methodology and discuss our results for the
multivariate normal case. In Appendix 6 we discuss further the generality of
this approach and show that it extends to other distributions of the
elliptical family including, in particular, the multivariate Student-t.

The logarithm of the likelihood for the normal case is proportional to:

$\ln\mathcal{L}(\bm{\theta;X}_{t})=\ln\left|\bm{J}\right|-(\bm{X}_{t}-\bm{\mu})\bm{J}(\bm{X}_{t}-\bm{\mu})^{T}+k\;\;,$
(3)

where $\bm{X}_{t}=(x_{t,1},x_{t,2},...,x_{t,n})$ is the $n$-dimensional
multivariate returns observation vector at time $t$; $\bm{\theta}$ is the
model parameter set, which includes $\bm{\mu}$ the vector of means and
$\bm{J}$ the generalized precision matrix and; $k$ is a constant which is
independent from ${\bm{\mu,J}}$ or $\bm{X}_{t}$ (see Appendix 6). In the
multivariate normal case $\bm{J}=\bm{\Sigma}^{-1}$ is the inverse of the
covariance. Our results generalize to other elliptical distributions with
defined covariance where $\bm{J}$ is proportional to the inverse covariance,
which we assume is defined and invertible. Results for the Student-t are
explicitly reported in appendix 6.

Our goal is to study the log-likelihood in Eq. 3 using different estimation
windows and comparing how the precision matrices, estimated through maximum-
likelihood and TMFG-LoGo, perform. We considered a dataset of daily closing
prices of US stocks entering among the constituents of the S&P 500 index
between 02/01/1997 and 31/12/2015. After screening for those continuously
traded and those not displaying abnormal returns, we reached a final dataset
of 342 stocks. For each asset $i=1,...,n$, we calculated the corresponding
daily returns $x_{t,i}=P_{t,i}/P_{(t-1),i}-1$, where $P_{t,i}$ is the closing
price of stock $i$ at time $t$, for a total of 4026 daily multivariate
observations.

We designed a resampling experiment in which we select 100 stocks at random
among the 342 and a random trading day spanned in our dataset. Starting from
the randomly selected trading day and going back in time, we define five train
sets of different sizes by including an increasing number of observations. We
start at 101 observations, then 150, 250, 500, 1000 and finally 1500
observations. We then use a fixed-length test set of 500 observations
following the randomly selected trading day. We keep the test set length fixed
to avoid biases and selected 500 observations so that, for all estimation
windows, the main crisis event (i.e. Global Financial Crisis in 2008) can be
randomly included in or out of sample. Figure 1 shows a sketched example of
our train/test split with different estimation windows.

Figure 1: Training and Testing scheme. We randomly sample the ending date of
the training period, the ‘trading day’. We then estimate the model parameters
considering different training windows using observations up to the randomly
selected trading day. The following 500 observations are used for testing.

We use the train set to estimate the mean vector $\bm{\mu}$, the maximum-
likelihood covariance matrix $\bm{\Sigma}$ and the sparse TMFG LoGo covariance
matrix $\bm{\Sigma_{TMFG}}$. These parameters are then used to compute the
log-likelihood in eq. 3 for both in-sample and out-of-sample observations. We
then investigate how the different estimates used in a portfolio optimization
procedure affect the optimal weights and portfolio characteristics. To this
extent, we considered the standard, unconstrained Markowitz optimization
problem described in Section 2.1. This is done to avoid any bias coming from
constraints in our analysis and to keep the framework as plain as possible. We
focus therefore our analysis on the minimum variance portfolio, that is the
efficient portfolio that minimizes the expected variance. To obtain the
solution for the minimum variance portfolio, the portfolio optimization
problem in eq. LABEL:min_variance rewrites:

$\displaystyle\min_{\bm{W}}$ $\displaystyle\sigma_{p}^{2}=\bm{W\Sigma W}^{T}$
(4) s.t. $\displaystyle\bm{W}\mathbbm{1}=1,$

which gives the optimal, minimum variance weights:

$\bm{W_{min}^{*}}=c\;\mathbbm{1}\bm{\Sigma^{-1}},$ (5)

where $c=\frac{1}{\mathbbm{1}^{T}\bm{\Sigma^{-1}}\mathbbm{1}}$ is a
normalization constant. Considering the estimation scheme described above and
outlined in fig. 1, the estimated covariance matrices are used as inputs in
the minimum variance optimal portfolio, eq. 5, to compare the different log-
likelihood levels obtained out-of-sample and the corresponding effects on
portfolio performances.

##  4 Results

###  4.1 Likelihood Comparison

Figure 2 reports the average log-likelihood for the train data (fig. 2(a)) and
for the test data (fig. 2(b)) computed across 100 resamplings. The larger the
log-likelihood is, the better the parameters $\bm{\theta}$ are at describing
the data, for the assumed model. Fig. 2(a) shows that, as expected and by
definition, the maximum-likelihood estimate of the covariance matrix provides
a higher in-sample likelihood as compared to the TMFG covariance, although the
latter tracks quite closely the maximum-likelihood. Also, one might observe
that the likelihood is strictly decreasing with the number of observations
included in the estimation window. Indeed, as the number of observations
decreases relative to the parameters, the model overfits the sample yielding
larger in-sample likelihoods. Filtering the covariance matrix and reducing the
number of parameters clearly limits the overfitting potential of the model as
shown by the lower levels of likelihood attained by TMFG when fewer
observations are used which, therefore, results in a larger gap in likelihood
relative to the maximum-likelihood covariance.

Perhaps more interestingly, fig. 2(b) reports the likelihoods obtained out-of-
sample using the two different in-sample estimates of the covariance matrix.
The first observation is that TMFG-LoGo provides a substantially larger log-
likelihood, especially for short estimation windows. This result is
exacerbated by the fact that when 101 observations are considered, the number
of stocks is very close to the number of observations in our samples. While
the resulting covariance is still full-rank (number of observations $>$ number
of variables), it leads to unstable estimation in the maximum-likelihood
covariance (i.e. the so-called “the curse of dimensionality”) whereas TMFG-
LoGo is still well defined. Note that there is a break y axis of the figure to
allow a better inspection of the results. The figure shows that for longer
estimation windows, the out-of-sample log-likelihood computed with the
maximum-likelihood covariance tends to converge to the TMFG likelihood which,
however, a) always provides the best out-of-sample likelihood in our
experiment and b) provides quite stable likelihood values also for shorter
estimation windows. We conclude that the TMFG-LoGo algorithm does a good job
at filtering the correlation structure providing higher out-of-sample
likelihood and stable results with shorter estimation windows, confirming the
results with stationary time series previously reported in (Barfuss et al.,
2016)

(a) Likelihood comparison in-sample

(b) Likelihood comparison out-of-sample. Note the y axis break to fit the
scale for 101 days estimation window.

Figure 2: Log Likelihood computed in sample and out of sample. Both the
likelihood computed using the maximum likelihood and the TMFG covariances
decrease in sample as the sample size increases. The maximum likelihood
covariance delivers by construction the highest likelihood, but the TMFG
likelihood tracks it closely. In test, instead, the TMFG covariance always
attain the highest likelihood and delivered good results also when the number
of observations becomes close to the number of variables.

###  4.2 Impact of precision matrix estimate on optimal portfolios

We now address empirically the question of what is the impact of different
parameter estimates on portfolios weights and performances, when these
parameters are used as inputs in the portfolio optimization problem in eq.
LABEL:min_variance. Having focused our attention on the minimum variance
portfolio on the efficient frontier, we report in Figure 3 the realized
standard deviation of portfolios obtained using the same parameters which
provided the log-likelihood displayed in Figure 2. The chart shows that,
overall, the out-of-sample portfolio variance decreases as the likelihood
increases up until when 750 observations are used. This is coherent with
respect to the likelihood results that reported, indeed, increasing
likelihoods for the same estimation windows. In particular, for shorter
estimation windows, the TMFG-LoGo covariance matrix provides portfolios with
significantly lower realized variance. Also, little changes are observed in
the realized variance when observations from 101 to 750 are included,
signalling that the TMFG-LoGo extract the relevant dependency links also when
few observations are available. The gap in performance tends to reduce as the
number of observations in the estimation window increases, with the TMFG-LoGo
portfolios always displaying lower volatility. However, when more than 750
daily observations are included, while the out of sample likelihood remains
flat or slightly increases, the portfolios’ variance tends to increase.

Figure 3: Realized Standard Deviation. Increasing the estimation window and
for higher values of Likelihood (figure 2), the realized standard deviation of
portfolios decreases. Y axis break to fit the scale for 101 days estimation
window..

To further investigate this pattern, we report in Figure 4 the volatilities
for all 100 resamplings and considering steps of 25 observations in the
estimation windows. The figure confirms that the TMFG-LoGo covariances
delivered overall less volatile portfolios across resamplings and estimation
windows. Secondly, the figure shows that the portfolios obtain the lowest out-
of-sample variance when approximately 2 to 3 years of daily observations (450
to 700 observations) are included in the train set. This pattern is clear for
the Maximum likelihood portfolios, with means, quintiles and outliers drifting
upwards when more than 750 observations are included. The TMFG filtered
covariance regularizes and smooths this effect as well, but still when more
than 750 observations are included, the resulting portfolios exhibit a
slightly higher variance. This is consistent with the literature showing that
longer estimation windows provide worse forecasts in financial time series due
to the regime-changing nature of financial markets (Procacci and Aste, 2019).
It is also worth emphasizing that the same mean vectors are used as inputs in
the Markowitz optimization for both the TMFG-LoGo and the Maximum-Likelihood
portfolios, hence the differences in performances are due solely to the
different estimates of the covariance matrix.

(a) Realized volatilities obtained with Maximum Likelihood covariances

(b) Realized volatilities obtained with TMFG filtered covariances

Figure 4: Portfolio realized volatility across resamplings for different
estimation windows. The box-plot shows the distribution of the variances
obtained for 100 resampled portfolios and the the blue line overlaid shows the
average variance (i.e. the mean of the variances’ distribution).

Finally, we address the impact of sparsity on portfolio construction by
analyzing what is the effect of a sparse precision matrix on the active bets
obtained from an optimization procedure. Figure 5 reports the number of Long
(fig. 6(a)) and Short (fig. 6(b)) positions (i.e. positive and negative
weights assigned to the stocks in portfolio) on average across the 100
resamplings. The first observation is that the number of long positions tends
to increase as the estimation window increases and coherently the short
positions diminish accordingly. Using TMFG-LoGo precision matrices anticipates
this behavior, in that these portfolios always display a greater number of
long positions also for short estimation windows. The intuition behind this
phenomenon is that over the long term, markets tend to post positive returns
and possible outliers in assets’ means and correlations are polished. This
intuition is confirmed by looking at the distribution of weights across
resamplings in Figure 6. This chart (note the different scales) shows that
using the TMFG-LoGo covariance matrix significantly improves the stability of
the optimal solutions, reducing outliers and avoiding “corner”, i.e. extreme
solutions which are a typical pitfall of the unconstrained Markowitz
optimization. This results shows that the correlation coefficient among assets
plays an important role in that for high correlation levels, the optimization
procedure would prefer one stock in place of another for slightly more
appealing mean or variance features. Having filtered the correlation structure
in the TMFG-LoGo procedure, we obtained a portfolio that is much more general
(hence the anticipated larger number of Long positions) and less sensitive to
single assets features given the filtered correlation among stocks.

(a) Number of “Buy” positions.

(b) Number of “Sell” positions.

Figure 5: Comparison of Buy/Sell Active Positions. As the number of training
observations increases, the optimizations delivers an increasing number “Long”
positions. This tendency is anticipated when using TMFG fiiltered covarinace
which always delivers an higher number of Long positions.

Lastly, considering the standard unconstrained optimization problem in eq.
LABEL:min_variance, both the maximum-likelihood and the TMFG matrices produce
portfolios that are in the vast majority of cases investing in all assets. In
other words, even considering a sparse precision matrix like in the TMFG-LoGo
case, we very rarely found weights equal to zero assigned to some assets.

(a) Distribution of optimal weights using Maximum-Likelihood covariance

(b) Distribution of optimal weights using TMFG covariance

Figure 6: Optimal Weights Distribution. Using the TMFG filtered covariance in
the optimization provides stable weights as compared to the maximum-likelihood
covariance, avoiding “corner solutions” and enhancing diversification.

###  4.3 Non Stationarity

From the results discussed in previous section and shown in Figure 3 and
Figure 2, we found that the portfolio performances improve coherently with the
likelihood up until approximately 3 years of observations are used in the
train set to estimate the models parameters. However, when more onservations
are included in the train set, the likelihood of the parameters detaches from
the portfolio performances and we speculate that this is due to the role of
nonstationarity. To further investigate this phenomenon, Figure 7 reports the
likelihood corresponding to each out-of-sample observation in our experiment.

(a) Mean Likelihood for each observation across resamplings. Comparison of
likelihoods obtained when 125, 750 and 1500 days are used in the train set.

(b) Boxplot of likelihoods representing the quartiles and min-max levels for
each observation across resamplings, having removed outliers. The plot is for
750 days train set (blue plot on the left).

Figure 7: Out-of-sample likelihood measured observation-by-observation.

Figure 7(a) shows the average likelihood across 100 resamplings for each out-
of-sample observation. We note that when shorter estimation windows are used
to estimate the models’ parameters (i.e. 125 days) the likelihood is higher in
the days immediately following the estimation window, but tends to rapidly
decrease as the observations depart from the training window. Larger
estimation windows (i.e. 750 or 1,500 days) instead, lead to a more stable
likelihood in the long run, but at the cost of a lower likelihood for the
obervations closer to the estimation set. Figure 7(b) shows the observation-
wise box-plot of the likelihood computed across the resamplings when 750
observation are used in the train set. The box plot reports the 25%-75%
quantile interval (dark blue) and the max-min interval (‘whiskers’ light blue)
having excluded the ‘outliers’ that are below the whiskers’ (Langford, 2006).
Other than decreasing means, the figure shows that as the observations depart
from the train set, the amount of observations posting a significantly lower
likelihood increases, together with the downside volatility. In other words,
it is more likely to have observations that are far from the model estimated
in sample, supporting the conclusions drawn from Figure 7(a).

As we discussed in (Procacci and Aste, 2019), market states tend to be
persistent in daily observations. Shorter estimation windows, therefore, tend
to better describe the system belonging to the same ‘state’ which is likely to
be persistent for adjacent observations. Notice that by considering the
aggregate behaviors across 100 resamplings, we want to avoid specific market
conditions and state shifts, but rather focus on the general behaviour. The
evolution of the financial system is obviously very dynamic and the goodness
of parameters is certainly dependent on both systemic and idiosyncratic
events.

These results also provide further insights on the findings discussed in
Figure 2 and Figure 3 in that our conclusions are dependent on the number of
out-of-sample observations that in our case coincides with the portfolio
holding period - i.e. 500 days in our experiments. Shorter estimation windows
provide better fit in the short term, while larger estimation windows provide
robustness in the long run. The optimal balance between these two effects
depends on the holding period and in our experiments it is achieved with
approximately 3 years observations in the estimation window. Short holding
periods do not require robustness in the long run (i.e. shorter estimation
windows would deliver better results). As the holding period increases, the
long term robustness becomes more relevant than the short term fit and larger
estimations windows have to be preferred.

##  5 Conclusions

Portfolio construction is a cornerstone of financial theory and practice.
However, it is still today a controversial topic for both academics and
practitioners. Any portfolio optimization strategy relies on assumptions and
modeling of the future market structure. However, inferring such structure
from past observations is a very challenging task, plagued by uncertainty
around parameters estimation and relying on some non fully satisfied
assumptions.

We identify three main sources of inaccuracies and errors: 1\. model
oversimplification; 2. limited size of the estimation set; 3. non-
stationarity. We address oversimplification by introducing a modeling that
uses a $L_{0}$-norm regularized elliptical multivariate distribution,
demonstrating that it over-performs traditional models both in likelihood and
in portfolio variance performances. We test the effect of sample size by
training the models on windows of different sizes and find that performances
initially increase with sample size but then eventually decrease for windows
above 750 days. We attribute the initial improvement in performance to
sampling error, which is reduced when more observations are included, and we
interpret the decay in performance when more the 750 observations are included
as an instance of non-stationarity. We further investigate this phenomenon by
studying the likelihood corresponding to individual observations out-of-sample
and show that shorter estimation windows deliver higher out-of-sample
likelihood in the days immediately following the train window, but it tends to
rapidly decrease afterwards. As more observations are included in the training
set, the out-of-sample likelihood gains stability, with larger values in the
long term, but at the cost of lower likelihood in the short term. We conclude
that the financial system changes significantly through time and the ‘optimal’
fit in finance needs to be defined in terms of the holding period.
Our main contribution to the literature on portfolio construction is the
demonstration of the relationship between the goodness of the model, measured
as out-of-sample likelihood, and the realized portfolio volatility. We show
that higher likelihood obtained with filtered TMFG-LoGo precision matrices
correspond to lower portfolio volatility out-of-sample. The relationship
between larger likelihood and lower realized volatility is also verified in
the maximum-likelihood estimate of the covariance matrix when computed over
train sets of different lengths. Further, we show that sparse, filtered
covariance matrices can significantly reduce estimation errors coming from
both sampling error and non-stationarity. It also reduces many of the
instability problems related to mean-variance optimal weights.

Finally, all of the analysis and conclusions drawn in this paper are based on
different estimates of the covariance matrix. While forecasting future returns
remains of primary importance in trading and wealth management, we showed that
the correlation structure, sometimes overlooked in the asset allocation
literature, plays a key role in portfolio construction and a good deal of
performances depend upon it.

## Acknowledgments

TA acknowledges partial support from ESRC (ES/K002309/1), EPSRC (EP/P031730/1)
and EC (H2020-ICT-2018-2 825215).

## References

  * Arditti and Levy (1976) Arditti, F. and Levy, H., Portfolio Efficiency Analysis in Three Moments: The Multiperiod Case. Journal of Finance, 1976, pp. 797–809.
  * Aste and Di Matteo (2006) Aste, T. and Di Matteo, T., Dynamical networks from correlations. , 2006, 370, 156–161.
  * Aste (2020) Aste, T., Topological regularization with information filtering networks. arXiv preprint arXiv:2005.04692, 2020.
  * Aste and Di Matteo (2017) Aste, T. and Di Matteo, T., Causality network retrieval from short time series. arXiv preprint arXiv:1706.01954, 2017.
  * Bamberg and Dorfleitner (2001) Bamberg, G. and Dorfleitner, G., Fat Tails and Traditional Capital Market Theory. Working Pape, University of Augsburg, 2001.
  * Barfuss et al. (2016) Barfuss, W., Massara, G.P., di Matteo, T. and Aste, T., Parsimonious modeling with Information Filtering Networks. Phys. Rev. E, 2016, 94, 062306.
  * Bawa (1978) Bawa, V.S., Safety-First, Stochastic Dominance, and Optimal Portfolio Choice. Journal of Financial and Quantitative Analysis, 1978, 13, 255–271.
  * Berkane and Bentler (1986) Berkane, M. and Bentler, P., Moments of elliptically distributed random variates. Statistics & Probability Letters, 1986, 4, 333 – 335.
  * Black and Litterman (1992) Black, F. and Litterman, R., Global Portfolio Optimization. Financial Analysts Journal, 1992, 48, 28–43.
  * Broadie (1993) Broadie, M., Computing efficient frontiers using estimated parameters. Annals of Operations Research, 1993, 45, 215–229.
  * Chamberlain (1983) Chamberlain, G., A characterization of the distributions that imply mean–Variance utility functions. Journal of Economic Theory, 1983, 29, 185–201.
  * Cramér (1946) Cramér, H., Mathematical Methods of Statistics, Princeton Landmarks in Mathematics and Physics, 1946, Princeton University Press.
  * Daniels (1965) Daniels, H.E., The asymptotic efficiency of a maximum likelihood estimator. Matematika, 1965, 9, 149–161.
  * De Franco et al. (2019) De Franco, C., Nicolle, J. and Pham, H., Bayesian Learning For The Markowitz Portfolio Selection Problem. International Journal of Theoretical and Applied Finance (IJTAF), 2019, 22, 1–40.
  * Fang et al. (1990) Fang, K., Kotz, S. and Ng, K., Symmetric multivariate and related distributions, no. 36 In Monographs on statistics and applied probability, 1990 (Chapman & Hall: London [u.a.]).
  * Friedman et al. (2008) Friedman, J., Hastie, T. and Tibshirani, R., Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 2008, 9, 432–441.
  * Gollier (2001) Gollier, C., Wealth Inequality and Asset Pricing. The Review of Economic Studies, 2001, 68, 181–203.
  * Harlow (1991) Harlow, W.V., Asset Allocation in a Downside-Risk Framework. Financial Analysts Journal, 1991, 47, 28–40.
  * Hogan and Warren (1974) Hogan, W.W. and Warren, J.M., Toward the Development of an Equilibrium Capital-Market Model Based on Semivariance. Journal of Financial and Quantitative Analysis, 1974, 9, 1–11.
  * Kritzman (2000) Kritzman, M.P., Puzzles of finance: six practical problems and their remarkable solutions, Wiley investment series, 2000 (John Wiley & Sons, Inc.: New York).
  * Langford (2006) Langford, E., Quartiles in Elementary Statistics. Journal of Statistics Education, 2006, 14.
  * Lauritzen (1996) Lauritzen, S.L., Graphical Models, 1996 (Oxford University Press: Oxford).
  * Ledoit and Wolf (2020) Ledoit, O. and Wolf, M., The Power of (Non-)Linear Shrinking: A Review and Guide to Covariance Matrix Estimation. Journal of Financial Econometrics, 2020.
  * Li and Ng (2000) Li, D. and Ng, W.L., Optimal Dynamic Portfolio Selection : Multi-period Mean-Variance Formulation. Mathematical Finance, 2000, 10, 387–406.
  * Lwin et al. (2017) Lwin, K.T., Qu, R. and MacCarthy, B.L., Mean-VaR portfolio optimization: A nonparametric approach. European Journal of Operational Research, 2017, 260, 751 – 766.
  * Mantegna (1999) Mantegna, R.N., Hierarchical structure in financial markets. The European Physical Journal B - Condensed Matter and Complex Systems, 1999, 11, 193–197.
  * Mantegna and Stanley (2000) Mantegna, R. and Stanley, H., An Introduction to Econophysics: Correlations and Complexity in Finance, 2000, Cambridge University Press.
  * Markowitz (1952) Markowitz, H., Portfolio Selection. The Journal of Finance, 1952, 7, 77–91.
  * Maruyama and Seo (2003) Maruyama, Y. and Seo, T., Estimation of Moment Parameter in Elliptical Distributions. Journal of the Japan Statistical Society. Japanese issue, 2003, 33, 215–229.
  * Massara and Aste (2019) Massara, G.P. and Aste, T., Learning clique forests. arXiv preprint arXiv:1905.02266, 2019.
  * Massara et al. (2015) Massara, G.P., di Matteo, T. and Aste, T., Network Filtering for Big Data: Triangulated Maximally Filtered Graph. CoRR, 2015, abs/1505.02445.
  * Michaud (1989) Michaud, R., The Markowitz Optimization Enigma: Is Optimized Optimal?. Working Paper, University of Augsburg, 1989, pp. 31–42.
  * Michaud and Michaud (1998) Michaud, R. and Michaud, R., Efficient Asset Management: A practical Guide to Stock Portfolio Optimization and Asset Allocation, 1998 (Harvard Business School Press: Boston).
  * Nawrocki (1996) Nawrocki, D., Portfolio Analysis with a Large Universe of Assets. Applied Economics, 1996, pp. 1191–1198.
  * Owen and Rabinovitch (1983) Owen, J. and Rabinovitch, R., On the Class of Elliptical Distributions and their Applications to the Theory of Portfolio Choice. The Journal of Finance, 1983, 38, 745–752.
  * Procacci and Aste (2019) Procacci, P.F. and Aste, T., Forecasting Market States. Quantitative Finance, 2019, 19, 1491–1498.
  * Scherer (2007) Scherer, B., Portfolio Construction and Risk Budgeting, 2007.
  * Scherer et al. (2012) Scherer, B., Winston, K. and O’Cinneide, C., Bayesian Methods In Investing, 2012, Oxford University Press.
  * Song et al. (2012) Song, W.M., Di Matteo, T. and Aste, T., Building complex networks with Platonic solids. Phys. Rev. E, 2012, 85, 046115.
  * Tumminello et al. (2005) Tumminello, M., Aste, T., Di Matteo, T. and Mantegna, R.N., A tool for filtering information in complex systems. Proceedings of the National Academy of Science, 2005, 102, 10421–10426.
  * Van Der Maaten et al. (2009) Van Der Maaten, L., Postma, E. and Van den Herik, J., Dimensionality reduction: a comparative review. Journal of Machine Learning Research, 2009, 10, 66–71.
  * Wald (1949) Wald, A., Note on the Consistency of the Maximum Likelihood Estimate. Annals of Mathematical Statistics, 1949, 20, 595–601.

##  6 Elliptical Distributions

Consider an $n$-dimensional vector of multivariate returns
$\bm{X}=(x_{1},x_{2},...,x_{n})$. If $\bm{X}$ is elliptical distributed, then
its probability density function is defined as:

$f_{X}(\bm{X})=c_{n}|\bm{J}|^{1/2}g_{n}\left[(\bm{X}-\bm{\mu})\bm{J}(\bm{X}-\bm{\mu})^{T}\right],$
(6)

where $\bm{\mu}\in\mathbb{R}^{1\times n}$ is the vector of location (mean)
parameters and $c_{n}$ is a normalization constant. The matrix,
$\bm{J}=\bm{\Omega}^{-1}\in\mathbb{R}^{n\times n}$ is the generalized
precision matrix, a positively defined matrix which is the inverse of the
dispersion matrix $\bm{\Omega}$. When the covariance is defined (as we assume
in this paper) then $\bm{\Omega}=(-\psi^{\prime}(0))^{-1}\bm{\Sigma}$, that
is, $\bm{\Omega}$ is proportional to the covariance matrix and the
proportionality factor is the inverse of the first derivative of the
characteristic generator evaluated at 0. The function, $g_{n}(\cdot)$ is
called density generator.

Also, let us stress that $(\bm{X}-\bm{\mu})\bm{J}(\bm{X}-\bm{\mu})^{T}$ \-
i.e. the generalized, square Mahalanobis distance - is a quadratic term and
hence a non-negative quantity provided that the matrix $\bm{\Omega}$ is
positive definite. To ease the notation, for the remaining of the paper we
shall refer to the generalized Mahalanobis distance as $d^{2}$:

$d^{2}=(\bm{X}-\bm{\mu})\bm{J}(\bm{X}-\bm{\mu})^{T}.$ (7)

For dfferent density generators $g_{n}(\cdot)$ we obtain different
distributions of the elliptical family. It is easy to see, for example, that
the normal distribution is obtained by using:

$g(u)=e^{-u/2},$ (8)

and $\bm{\Omega}=\bm{\Sigma}$.

Similarly the Student-t distribution is obtained by using:

$g_{n}(u)=\left(1+\frac{u}{v}\right)^{-\frac{n+v}{2}},$ (9)

where $v$ is the degrees of freedom, and
$\bm{\Omega}=\frac{\nu}{\nu-2}\bm{\Sigma}$.

The validity of the mean-variance framework for elliptical distributions has
long been established in literature (Owen and Rabinovitch, 1983). This
proposition is derived easily from two properties of the elliptical
distributions. First, for every elliptical distribution with defined mean and
variance, the distribution is completely specified by them (Owen and
Rabinovitch (1983) or Chamberlain (1983)), with all the higher moments being
either zero or proportional to the first or second moment. Second, any linear
combination of multivariate elliptically distributed variables is also an
elliptically distributed variable. In the case of normal distribution and
Stiudent-t distribution they also have the same density generator function.
Further details on these properties are provided in Appendix 7.

It follows that, if asset returns have a multivariate elliptical distribution
$\bm{X}\sim\mathcal{E}_{n}(\bm{\mu,\Omega},g_{n})$, then the portfolio
expected return and dispersion are given by, respectively,
$\mathbb{E}[r_{p}]=\bm{W\mu}^{T}$ and $\sigma_{p}=\bm{W\Omega W^{T}}$,
matching the optimization framework outlined in Section 2.1.

With respect to our likelihood analysis, considering distributions with
probability density function of the form specified in Eq. 6, the corresponding
likelihood function is of the form:

$\mathcal{L}_{ED}(\bm{\theta;X})=\left|\bm{J}\right|^{1/2}\;g_{n}\left(d^{2}\right)\;.$
(10)

Where $ED$ denotes the general Elliptical Distributions and we omitted the
constant of integration.

To stress the general validity of our analysis for other elliptical
distributions, we repeated the experiments discussed in Section 3 considering
the t-student generator.

Assuming a Student - t distribution of the log returns, the log likelihood
(Eq. 10) is:

$\ln\mathcal{L}_{Student}=\frac{\ln|\bm{J}|}{2}-\frac{n+\nu}{2}\ln\left(1+\frac{d^{2}}{\nu-2}\right)$
(11)

where $n$ is the sample size and $\nu$ is the degree of freedom. Figure 8(a)
reports the likelihood comparison for the same resamplings as in Figure 2 but
using a student-t log likelihood as in Eq. 11. Here we used $n=500$
observations (i.e. the out-of-sample size) and $\nu=3$. We verified that this
findings are robust across different degrees of freedom in the range
$\nu=[2.1,4]$.

(a) Log likelihood assuming a t-student distribution of log returns. Figure 8:
Out-of-sample log likelihood likelihood computed using the maximum likelihood
and the TMFG covariances. These results are coherent with the findings related
to the Normal distribution presented in Section 4.1

##  7 Properties of Elliptical Distributions

In this section we recall some useful properties of Elliptical Distribution
which we referred to in our discussion and particularly in Section 6.

######  Property 1 (Distribution Definition).

Consider an n-dimensional random vector $\bm{X}=(X_{1},...,X_{n})$. $\bm{X}$
has a multivariate elliptical distribution with location parameter $\bm{\mu}$
and dispersion parameter $\bf\Omega$,written as
$\bm{X}\sim\mathcal{E}(\bm{\mu,\Omega})$ if its characteristic function $\phi$
can be expressed as:

$\phi_{X}(\bm{w})=\mathbb{E}(e^{i\bm{w}\bm{X}})=e^{i\bm{w}\bm{\mu}}\psi\left({\frac{1}{2}\bm{w}\bm{\Omega}\bm{w}^{T}}\right),$
(12)

for some location parameter $\bm{\mu}\in\mathbb{R}^{1\times n}$, positive-
definite dispersion matrix $\bm{\Omega}\in\mathbb{R}^{n\times n}$ and for some
function $\psi(\cdot):[0,\infty)\rightarrow\mathbb{R}$ such that
$\psi\left(\sum_{i=1}^{n}w^{2}_{i}\right)$ is a characteristic function, which
is called characteristic generator. If
$\bm{X}\sim\mathcal{E}(\bm{\mu,\Omega})$ and if its density $f_{X}(\bm{X})$
exists, it is of the form defined in Eq. 6.

######  Property 2 (Density Generator).

The function $g(\cdot)$ defined in Section 6 is guaranteed to be density
generator if the following condition holds:

$\int_{0}^{\infty}x^{n/2-1}g_{n}(x)dx<\infty.$ (13)

######  Property 3 (Affine Equivariance).

If $\bm{X}=(X_{1},...,X_{n})$ is an n-dimensional elliptical random variable
with location parameter $\bm{\mu}$ and dispersion parameter $\bm{\Omega}$ so
that $\bm{X}\sim\mathcal{E}_{X}(\bm{\mu,\Omega})$, then for any vector
$\bm{a}\in\mathbb{R}^{1\times m}$ and any matrix $\bm{B}\in\mathbb{R}^{m\times
n}$ the following affine equivariace holds:

$\bm{Y}=\bm{a+BX}\sim\mathcal{E}_{Y}(\bm{a+B\mu,B\Omega B}).$ (14)

In other words, any linear combination of multivariate elliptical
distributions is another elliptical distribution. In the special cases of
normal, Student-t and Cauchy distributions, the induced density generators are
$m$-dimensional version of the original generator of $\bm{X}$.

For the proof of Properties 1,2 and 3, we refer to Fang et al. (1990).

This implies that any portfolio $Y=\beta_{1}X_{1}+...+\beta_{n}X_{n}$ of
elliptically distributed variables is distributed accordingly with a
(univariate) elliptical distribution, which is a location-scale distribution.
Furthermore, for any univariate elliptical distribution all moments can be
obtained from the first and second moments (if defined). In particular, for
centered variables with zero mean ($\mu_{Y}=0$), the resulting distribution of
$Y$ is symmetrical around zero and it has all odd moments equal to zero and
all even moments given by:

$\mu_{2m}=c_{m}\mu_{2}^{m},$

with

$c_{m}=\frac{(2m)!}{(2^{m}m!)}\frac{\psi^{(m)}(0)}{(\psi^{(1)}(0))^{m}}.$

Where $\psi^{(m)}(0)$ indicated the $m^{th}$ derivative of $\psi(\omega)$
computed at $\omega=0$.

As an example, in the normal (0,1) case, $\mu_{2}=1$, $c_{m}=0$ for all
$m=1,2,...$, the kurtosis is $\mu_{(4)}=\frac{4!}{2^{4}4!}=3$, and
$\mu_{(2m)}=\frac{(2m)!}{(2^{m}m!)}$. For the proof we refer to Berkane and
Bentler (1986), which derived this property by succesive differentiations of
$\phi(\cdot)$, and to Maruyama and Seo (2003), which attained the same result
by expressing the elliptical distribution in terms of a random vector with
uniform distribution on the unit sphere.

Therefore the mean-variance optimization is of general applicability and
relevance for any portfolio generated from multivariate elliptically
distributed variables.


$\int_{Q}(a_{1}\overline{u_{1}})\partial_{t}udV_{g}dt=-\int_{Q}\partial_{t}(a_{1}\overline{u_{1}})udV_{g}dt=-\int_{Q}(\partial_{t}a_{1}\overline{u_{1}}+a_{1}\partial_{t}\overline{u_{1}})udV_{g}dt.$
(5.4)

In particular, there are no boundary terms since $u|_{t=0}=u|_{t=T}=0$.

We then multiply equation (5.1) by $\overline{u_{1}}$ and integrate over $Q$.
Therefore, we deduce from (5.1)–(5.4) that

$\displaystyle\int_{Q}(a\partial_{t}u_{2}+qu_{2})\overline{u_{1}}dV_{g}dt$
(5.5)
$\displaystyle=\langle\mathcal{L}_{c,g,a_{1},q_{1}}u,\overline{u_{1}}\rangle_{L^{2}(Q)}-\langle
u,\overline{\mathcal{L}_{c,g,a_{1},q_{1}}^{*}u_{1}}\rangle_{L^{2}(Q)}$
$\displaystyle=\langle(c^{-1}\partial_{t}u,-\nabla_{g}u)\cdot\bar{\nu},\overline{u_{1}}\rangle_{L^{2}(\partial
Q)}-\langle
u,(c^{-1}\partial_{t}\overline{u_{1}},-\nabla_{g}\overline{u_{1}})\cdot\bar{\nu}\rangle_{H^{1/2}(\partial
Q),H^{-1/2}(\partial Q)}.$

Since $u|_{\Sigma}=u|_{t=0}=u|_{t=T}=0$, the second term on the right-hand
side of (5.5) vanishes. Furthermore, since $\partial_{\nu}u\in L^{2}(\Sigma)$
with $\partial_{t}u|_{t=0}=\partial_{\nu}u|_{V}=0$, we obtain the integral
identity

$\displaystyle\int_{Q}(a\partial_{t}u_{2}+qu_{2})\overline{u_{1}}dV_{g}dt=-\int_{\Sigma\setminus
V}\partial_{\nu}u\overline{u_{1}}dS_{g}dt+\int_{M}c^{-1}\partial_{t}u(T,x)\overline{u_{1}(T,x)}dV_{g}.$
(5.6)

We shall next substitute the CGO solutions (4.77) and (4.78) into (5.6),
multiply the equation by $h$, and pass to the limit $h\to 0$. In order to
analyze the limit of the terms on the left-hand side of (5.6), we use
estimates (3.43) and (4.4) to obtain the following estimates for the remainder
terms:

$\|r_{j}\|_{L^{2}(Q)}\leq\|r_{j}\|_{H^{1}_{\mathrm{scl}}(Q^{\mathrm{int}})}=o(1),\quad
j=1,2,$ (5.7)

and

$\|\partial_{t}r_{j}\|_{L^{2}(Q)}\leq\frac{1}{h}\|r_{j}\|_{H^{1}_{\mathrm{scl}}(Q^{\mathrm{int}})}=o(h^{-1}),\quad
j=1,2.$ (5.8)

On the other hand, the following lemma explains the behaviors of the two terms
on the right-hand side of (5.6) as $h\to 0$.

######  Lemma 5.1.

Let $u_{1}$ and $u$ be the functions described above. Then the following
estimates hold as $h\to 0$:

$\int_{M}c^{-1}\partial_{t}u(T,x)\overline{u_{1}(T,x)}dV_{g}=\left\\{\begin{array}[]{cc}\mathcal{O}(h^{-1/2}),&\text{if
}a\neq 0,\\\ \\\ \mathcal{O}(h^{1/2}),&\text{ if }a=0.\end{array}\right.$
(5.9) $\int_{\Sigma\setminus
V}\partial_{\nu}u\overline{u_{1}}dS_{g}dt=\left\\{\begin{array}[]{cc}o(h^{-1}),&\text{if
}a\neq 0,\\\ \\\ o(1),&\text{ if }a=0.\end{array}\right.$ (5.10)

We will postpone the proof of this result and use it to prove the uniqueness
of the damping coefficient first.

###  5.1. Uniqueness of the damping coefficient

From the respective CGO solutions (4.77) and (4.78) for $u_{1}$ and $u_{2}$,
we compute that

$u_{2}\overline{u_{1}}=e^{2i\lambda(\beta
t+x_{1})}c^{-\frac{n-2}{2}}(\overline{v_{s}}w_{s}+\overline{v_{s}}r_{2}+\overline{r_{1}}w_{s}+\overline{r_{1}}r_{2}).$

Therefore, by estimates (4.4), (4.5), (5.7), and the Cauchy-Schwartz
inequality, we have

$\displaystyle\|u_{2}\overline{u_{1}}\|_{L^{1}(Q)}\leq\mathcal{O}(1).$

Hence, the following estimate holds

$h\bigg{|}\int_{Q}qu_{2}\overline{u_{1}}dV_{g}dt\bigg{|}\leq
h\|q\|_{L^{\infty}(Q)}\|u_{2}\overline{u_{1}}\|_{L^{1}(Q)}=\mathcal{O}(h),\quad
h\to 0.$ (5.11)

We next consider the term
$h\int_{Q}a(\partial_{t}u_{2})\overline{u_{1}}dV_{g}dt$ on the left-hand side
of (5.6). To that end, direct computations yield

$\displaystyle(\partial_{t}u_{2})\overline{u_{1}}=$ $\displaystyle
e^{2i\lambda(\beta
t+x_{1})}c^{-\frac{n-2}{2}}[(\overline{v_{s}}\partial_{t}w_{s}+\overline{v_{s}}\partial_{t}r_{2}+\overline{r_{1}}\partial_{t}w_{s}+\overline{r_{1}}\partial_{t}r_{2})$
$\displaystyle+s\beta(\overline{v_{s}}w_{s}+\overline{v_{s}}r_{2}+\overline{r_{1}}w_{s}+\overline{r_{1}}r_{2})].$

Using estimates (4.4), (4.5), (5.7), (5.8), as well as the Cauchy-Schwartz
inequality, we obtain

$h\bigg{|}\int_{Q}e^{2i\lambda(\beta
t+x_{1})}c^{-\frac{n-2}{2}}a(\overline{v_{s}}\partial_{t}w_{s}+\overline{v_{s}}\partial_{t}r_{2}+\overline{r_{1}}\partial_{t}w_{s}+\overline{r_{1}}\partial_{t}r_{2})dV_{g}dt\bigg{|}=o(1),\quad
h\to 0,$

and

$h\bigg{|}\int_{Q}e^{2i\lambda(\beta
t+x_{1})}c^{-\frac{n-2}{2}}as\beta(\overline{v_{s}}r_{2}+\overline{r_{1}}w_{s}+\overline{r_{1}}r_{2})dV_{g}dt\bigg{|}=o(1),\quad
h\to 0.$

Therefore, we have

$h\int_{Q}a\partial_{t}u_{2}\overline{u_{1}}dV_{g}dt\to\int_{Q}e^{2i\lambda(\beta
t+x_{1})}c^{-\frac{n-2}{2}}\beta a\overline{v_{s}}w_{s}dV_{g}dt,\quad h\to 0.$
(5.12)

Using (5.11), (5.12), and Lemma 5.1, we deduce from (5.6) that

$\int_{Q}e^{2i\lambda(\beta t+x_{1})}c^{-\frac{n-2}{2}}\beta
a\overline{v_{s}}w_{s}dV_{g}dt\to 0,\quad h\to 0.$

On the other hand, since $a_{1},a_{2}\in W^{1,\infty}(Q)$ and $a_{1}=a_{2}$ on
the boundary $\partial Q$, we can continuously extend $a$ on
$(\mathbb{R}^{2}\times M_{0})\setminus Q$ by $0$ and denote the extension by
the same letter. Using $dV_{g}=c^{\frac{n}{2}}dV_{g_{0}}dx_{1}$, the change of
coordinates (4.26), Fubini’s theorem, the dominated convergence theorem, and
the concentration property (4.57) of the quasimodes $v_{s}$ and $w_{s}$, we
obtain

$\displaystyle\int_{Q}e^{2i\lambda(\beta t+x_{1})}c^{-\frac{n-2}{2}}\beta
a\overline{v_{s}}w_{s}dV_{g}dt$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}}\int_{\mathbb{R}}\int_{M_{0}}e^{2i\lambda(\beta
t+x_{1})}c\beta a\overline{v_{s}}w_{s}dV_{g_{0}}dx_{1}dt$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}}\int_{\mathbb{R}}\int_{M_{0}}\beta^{2}e^{2i\lambda((\beta^{2}-1)\widetilde{t}+p)}(ca)\overline{v_{s}}w_{s}dV_{g_{0}}dpd\widetilde{t}$
$\displaystyle\to$
$\displaystyle\beta^{2}(1-\beta^{2})^{-\frac{n-6}{4}}\int_{\mathbb{R}}\int_{\mathbb{R}}\int_{0}^{\frac{L}{\sqrt{1-\beta^{2}}}}e^{2i\lambda((\beta^{2}-1)\widetilde{t}+p)-2(1-\beta^{2})\lambda
r}(ca)(\beta\widetilde{t},p-\widetilde{t},\gamma(\sqrt{1-\beta^{2}}r))$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\times
e^{\overline{\Phi_{1}(\widetilde{t},p,r)}+\Phi_{2}(\widetilde{t},p,r)}\eta(\widetilde{t},p,r)drdpd\widetilde{t},\quad
h\to 0.$

As $\beta\in[\frac{1}{\sqrt{3}},1)$, we conclude that

$\displaystyle\int_{0}^{\frac{L}{\sqrt{1-\beta^{2}}}}\int_{\mathbb{R}}\int_{\mathbb{R}}$
$\displaystyle
e^{2i\lambda((\beta^{2}-1)\widetilde{t}+p)-2(1-\beta^{2})\lambda
r}(ca)(\beta\widetilde{t},p-\widetilde{t},\gamma(\sqrt{1-\beta^{2}}r))$ (5.13)
$\displaystyle\times
e^{\overline{\Phi_{1}(\widetilde{t},p,r)}+\Phi_{2}(\widetilde{t},p,r)}\eta(\widetilde{t},p,r)dpd\widetilde{t}dr=0.$

We next follow the arguments in [46, Section 4] closely to prove that (5.13)
holds when
$e^{\overline{\Phi_{1}(\widetilde{t},p,r)}+\Phi_{2}(\widetilde{t},p,r)}\eta(\widetilde{t},p,r)$
is removed from the integral. To this end, let us write
$\eta(\widetilde{t},p,r):=\eta_{1}(\tilde{t},r)\eta_{2}(p)$ in (4.58) with
$\partial\eta_{1}=0$, where
$\partial=\frac{1}{2}(\partial_{\widetilde{t}}+i\partial_{r})$, and denote

$\Psi(\widetilde{t},p,r)=e^{2i\lambda[(\beta^{2}-1)\widetilde{t}+p-i(\beta^{2}-1)r]}\eta_{1}(\widetilde{t},r).$
(5.14)

It follows from direct computations that
$\partial\Psi=\frac{1}{2}(\partial_{\widetilde{t}}-i\partial_{r})\Psi=0$.

Since $a$ is supported in $Q$, we get from (5.13) that

$\int_{\mathbb{R}}\int_{\mathbb{R}}\int_{\mathbb{R}}\Psi(\widetilde{t},p,r)(ca)(\beta\widetilde{t},p-\widetilde{t},\gamma(\sqrt{1-\beta^{2}}r))e^{\overline{\Phi_{1}(\widetilde{t},p,r)}+\Phi_{2}(\widetilde{t},p,r)}\eta_{2}(p)dpd\widetilde{t}dr=0,$

and there exists a constant $R>0$ such that $\hbox{supp }a\subset\subset
B_{\widetilde{t},p,r}(0,R)$. Also, since $\eta_{2}\in C^{\infty}(\mathbb{R})$
is arbitrary, for almost every $p\in\mathbb{R}$ we have

$\displaystyle\int_{\Omega_{p}}\Psi(\widetilde{t},p,r)(ca)(\beta\widetilde{t},p-\widetilde{t},\gamma(\sqrt{1-\beta^{2}}r))e^{\overline{\Phi_{1}(\widetilde{t},p,r)}+\Phi_{2}(\widetilde{t},p,r)}d\widetilde{t}dr=0,$
(5.15)

where $\Omega_{p}=\\{(\widetilde{t},r):(\widetilde{t},p,r,y)\in Q\\}$. We
shall view $\Omega_{p}$ as a domain in the complex plane with the complex
variable $z=\widetilde{t}+ir$.

Recall that it was explained in the formulae (3.1)–(3.4) how to transform the
hyperbolic operator $\mathcal{L}_{c,g,a,q}$ into another operator
$\mathcal{L}_{\widetilde{g},\widetilde{a},\widetilde{q}}$ of the same type,
where the contribution of the conformal factor $c$ was moved from the highest
order term to the lower order ones. In particular, we have $\widetilde{a}=ca$,
and the construction of the Gaussian beams in Section 4 is carried over for
this damping coefficient. Hence, equation (4.55) yields

$\partial(\overline{\Phi_{1}}+\Phi_{2})=\frac{1}{4}\beta(ca)(\beta\widetilde{t},p-\widetilde{t},\gamma(\sqrt{1-\beta^{2}}r)),$
(5.16)

Thus, it follows from (5.15) and (5.16) that

$\int_{\Omega_{p}}\partial\big{(}\Psi(z,p)e^{\overline{\Phi_{1}(z,p)}+\Phi_{2}(z,p)}\big{)}dz\wedge
d\overline{z}=0\quad\text{for almost every }p.$ (5.17)

We now discuss the regularity of $\Phi_{i}$. To that end, using (5.16) along
with the fact that $a(\cdot,p)\in L^{\infty}(\mathbb{C})$ is compactly
supported for almost every $p$, we see that $\partial\Phi_{i}\in
L^{p}(\mathbb{C})$ for $1\leq p\leq\infty$. By the boundedness of the
Beurling-Ahlfors operator $\overline{\partial}\partial^{-1}$ on
$L^{p}(\mathbb{C})$, $1<p<\infty$, we get
$\overline{\partial}\Phi_{i}=\overline{\partial}\partial^{-1}(\partial\Phi_{i})$,
which implies that $\nabla_{g}\Phi_{i}\in L^{p}(\mathbb{C})$. Furthermore,
since $\Phi_{i}\in L^{\infty}(\mathbb{C})$, we have $\Phi_{i}(\cdot,p)\in
W^{1,p}_{\mathrm{loc}}(\mathbb{C})$, $1<p<\infty$. Hence, we conclude that
$\Phi_{i}(\cdot,p)\in H^{1}(\Omega_{p})$, $i=1,2$. Thus, an application of
Stokes’ theorem [67, Theorem 18A] yields

$\int_{\partial\Omega_{p}}\Psi(z,p)e^{\overline{\Phi_{1}(z,p)}+\Phi_{2}(z,p)}d\overline{z}=0.$
(5.18)

By [24, Lemma 5.1], see also [2, Lemma 3, Section 2.3, Chapter 4], there
exists a non-vanishing function $F\in C(\overline{\Omega_{p}})$, anti-
holomorphic in $\Omega_{p}$, such that

$F|_{\partial\Omega_{p}}=e^{\overline{\Phi_{1}}+\Phi_{2}}|_{\partial\Omega_{p}}.$

Furthermore, the arguments in the proof of [43, Lemma 7.3] show that there
exists an anti-holomorphic function $G\in C(\overline{\Omega_{p}})$ such that
$F=e^{G}$ in $\Omega_{p}$, and we may assume that
$G=\overline{\Phi_{1}}+\Phi_{2}$ on $\partial\Omega_{p}$. Choosing
$\eta_{1}=Ge^{-G}$ in (5.14), we get from (5.18) that

$\int_{\partial\Omega_{p}}\big{(}\overline{\Phi_{1}(\widetilde{t},p,r)}+\Phi_{2}(\widetilde{t},p,r)\big{)}e^{2i\lambda[(\beta^{2}-1)\widetilde{t}+p-i(\beta^{2}-1)r]}d\overline{z}=0.$

Applying Stokes’ theorem again, using (5.16), and integrating over the $p$
variable, we obtain

$\int_{0}^{\frac{L}{\sqrt{1-\beta^{2}}}}\int_{\mathbb{R}}\int_{\mathbb{R}}e^{2i\lambda((\beta^{2}-1)\widetilde{t}+p)-2(1-\beta^{2})\lambda
r}(ca)(\beta\widetilde{t},p-\widetilde{t},\gamma(\sqrt{1-\beta^{2}}r))dpd\widetilde{t}dr=0.$
(5.19)

Finally, we use (4.26) to return to $(t,x_{1},\tau)$ coordinates from
$(\widetilde{t},p,r)$ coordinates and replace $2\lambda$ with $\lambda$. After
these changes, (5.19) becomes

$\int_{0}^{L}\int_{\mathbb{R}}\int_{\mathbb{R}}e^{i\lambda(\beta
t+x_{1})-\sqrt{1-\beta^{2}}\lambda\tau}(ca)(t,x_{1},\gamma(\tau))dx_{1}dtd\tau=0.$
(5.20)

We are now ready to utilize Assumption 1.1, the invertibility of the
attenuated geodesic ray transform on $(M_{0},g_{0})$. To that end, we let
$\mathcal{F}_{(t,x_{1})\to(\xi_{1},\xi_{2})}$ be the Fourier transform in the
Euclidean variables $(t,x_{1})$ and define

$\displaystyle f(x^{\prime},\beta,\lambda)$
$\displaystyle:=\int_{\mathbb{R}}\int_{\mathbb{R}}e^{i\lambda(\beta
t+x_{1})}(ca)(t,x_{1},x^{\prime})dx_{1}dt$
$\displaystyle=\mathcal{F}_{(t,x_{1})\to(\xi_{1},\xi_{2})}(ca)|_{(\xi_{1},\xi_{2})=-\lambda(\beta,1)}\quad\text{for
}x^{\prime}\in M_{0},\>\beta\in[1/2,1),\>\lambda\in\mathbb{R}.$

Since $a\in W^{1,\infty}(Q)$, we see that the function
$f(\cdot,\beta,\lambda)$ is continuous on $M_{0}$. Furthermore, as $\gamma$ is
an arbitrarily chosen non-tangential geodesic in $(M_{0},g_{0})$, we get from
(5.20) that the following attenuated geodesic ray transform vanishes,

$\int_{0}^{L}e^{-\sqrt{1-\beta^{2}}\lambda\tau}f(\gamma(\tau),\beta,\lambda)d\tau=0.$
(5.21)

By Assumption 1.1, there exists $\varepsilon>0$ such that
$f(\gamma(\tau),\beta,\lambda)=0$ whenever
$\sqrt{1-\beta^{2}}|\lambda|<\varepsilon$. Hence, there exist
$\beta_{0}\in(\frac{1}{\sqrt{3}},1)$, $\lambda_{0}>0$, and $\delta>0$ such
that for every $(\lambda,\beta)\in\mathbb{R}^{2}$ that satisfies
$|\beta-\beta_{0}|$, $|\lambda-\lambda_{0}|<\delta$, and $\lambda\neq 0$, we
have $\sqrt{1-\beta^{2}}|\lambda|<\varepsilon$. In particular, the mapping
$(\lambda,\beta)\mapsto-\lambda(\beta,1)$ is a diffeomorphism when
$\lambda\neq 0$, implying that
$\mathcal{F}_{(t,x_{1})\to(\xi_{1},\xi_{2})}(ca)=0$ in an open set of
$\mathbb{R}^{2}$. Lastly, the compact support of $a$ and Paley-Wiener theorem
yield that $\mathcal{F}(ca)$ is real analytic. Therefore, we get $ca=0$ in
$Q$. Since $c$ is a positive function, we must have $a=a_{1}-a_{2}=0$.

To complete the proof of uniqueness for the damping coefficient in Theorem
1.2, we still need to verify Lemma 5.1.

###### Proof of Lemma 5.1.

Let us prove estimate (5.9) first. To this end, using estimates (4.4) and
(5.7), the CGO solution (4.77), and the Cauchy-Schwartz inequality, we get

$\displaystyle\bigg{|}\int_{M}c^{-1}\partial_{t}u(T,x)\overline{u_{1}(T,x)}dV_{g}\bigg{|}$
$\displaystyle\leq\|c^{-1}\|_{L^{\infty}(Q)}\int_{M}|\partial_{t}u(T,x)\overline{u_{1}(T,x)}|dV_{g}$
(5.22) $\displaystyle\leq\mathcal{O}(1)\int_{M}|c^{-\frac{n-2}{4}}e^{-s(\beta
T+x_{1})}\partial_{t}u(T,x)|(|v_{s}(T,x)|+|r_{1}(T,x)|)dV_{g}$
$\displaystyle\leq\mathcal{O}(1)\|e^{-s(\beta
T+x_{1})}\partial_{t}u(T,\cdot)\|_{L^{2}(M)}.$

Utilizing the boundary Carleman estimate (3.23) and equation (5.1), we obtain

$\displaystyle\|e^{-s(\beta
T+x_{1})}\partial_{t}u(T,\cdot)\|_{L^{2}(M)}\leq\mathcal{O}(h^{1/2})\|e^{-s(\beta
t+x_{1})}\mathcal{L}_{c,g,a_{1},q_{1}}u\|_{L^{2}(Q)}.$ (5.23)

We then substitute the CGO solution (4.78) for $u_{2}$ into (5.1) to get

$\displaystyle e^{-s(\beta t+x_{1})}\mathcal{L}_{c,g,a_{1},q_{1}}u$
$\displaystyle=e^{-s(\beta t+x_{1})}(a\partial_{t}u_{2}+qu_{2})$
$\displaystyle=c^{-\frac{n-2}{2}}[as\beta(w_{s}+r_{2})+a(\partial_{t}w_{s}+\partial_{t}r_{2})+q(w_{s}+r_{2})].$

We recall that $s=h^{-1}+i\lambda$ and use estimates (4.5), (5.7), and (5.8)
to obtain

$\|e^{-s(\beta
t+x_{1})}\mathcal{L}_{c,g,a_{1},q_{1}}u\|_{L^{2}(Q)}=\left\\{\begin{array}[]{cc}\mathcal{O}(h^{-1}),&\text{if
}a\neq 0,\\\ \\\ \mathcal{O}(1),&\text{ if }a=0.\end{array}\right.$ (5.24)

Therefore, estimates (5.22)–(5.24) imply (5.9).

We next prove estimate (5.10). To that end, for all $\varepsilon>0$ we set

$\partial M_{+,\varepsilon}=\\{x\in\partial
M:\partial_{\nu}\varphi(x)>\varepsilon\\},\quad\text{and}\quad\Sigma_{+,\varepsilon}=(0,T)\times\partial
M_{+,\varepsilon}.$

We recall that in Section 1 we defined the open sets
$U^{\prime},V^{\prime}\subset\partial M$ such that they contain the back and
front faces $\partial M_{+},\partial M_{-}$ of the manifold $M$, respectively.
By the compactness of $\\{x\in\partial M:\partial_{\nu}\varphi(x)=0\\}$, there
exists $\varepsilon>0$ such that $\Sigma\setminus
V\subset\Sigma_{+,\varepsilon}$, where we had set $V=(0,T)\times V^{\prime}$.

We utilize estimates (4.4), (4.77), (5.7), as well as the Cauchy-Schwartz
inequality to get

$\displaystyle\bigg{|}\int_{\Sigma\setminus
V}\partial_{\nu}u\overline{u_{1}}dS_{g}dt\bigg{|}$
$\displaystyle\leq\int_{\Sigma_{+,\varepsilon}}e^{-s(\beta
t+x_{1})}|\partial_{\nu}u|(|v_{s}|+|r_{1}|)dS_{g}dt$ $\displaystyle\leq
C\bigg{(}\int_{\Sigma_{+,\varepsilon}}|e^{-s(\beta
t+x_{1})}\partial_{\nu}u|^{2}dS_{g}dt\bigg{)}^{1/2}(\|v_{s}\|_{L^{2}(\Sigma_{+,\varepsilon})}+\|r_{1}\|_{L^{2}(\Sigma)}).$

Next we estimate the terms in the inequality above. By Proposition 4.5 and
estimate (5.7), in conjunction with the inequalities

$\|r_{1}\|_{L^{2}(\Sigma)}\leq\|r_{1}\|_{L^{2}(Q)}^{1/2}\|r_{1}\|_{H^{1}(Q)}^{1/2}\quad\text{and}\quad\|r_{1}\|_{H^{1}(Q)}\leq
Ch^{-1}\|r_{1}\|_{H^{1}_{\mathrm{scl}}(Q)},$

we obtain

$\|r_{1}\|_{L^{2}(\Sigma)}=o(h^{-1/2}),\quad h\to 0.$ (5.25)

We now follow the steps in the proof of [21, Theorem 6.2] to verify

$\|v_{s}\|_{L^{2}(\Sigma_{+,\varepsilon})}=\mathcal{O}(1),\quad h\to 0.$
(5.26)

Due to the product structure of $\Sigma_{+,\varepsilon}$ and the fact that
$T>0$ is finite, it suffices to prove that $\|v_{s}\|_{L^{2}(\partial
M_{+,\varepsilon})}=\mathcal{O}(1)$.

Let $\pi_{1}:\partial M\to\mathbb{R}$ be a projection defined by

$\pi_{1}(x)=x_{1},\quad\text{for any }x=(x_{1},x^{\prime})\in\partial M.$

Without loss of generality, we assume that $\partial M_{\varepsilon}:=\partial
M\cap\pi^{-1}_{1}(\varepsilon)$ is a manifold for any $\varepsilon>0$. Hence,
$\partial M_{+,\varepsilon}$ is a compact $(n-1)$-dimensional manifold with
boundary.

First, we observe that $\partial M_{+,\varepsilon}\subset\partial M$ is an
open and precompact manifold of dimension $n-1$, same as $M_{0}$. Clearly, the
projections $\pi_{1}\colon\partial
M\to\mathbb{R},\>\pi_{1}(x_{1},x^{\prime})=x_{1}$, and $\pi_{2}:\partial M\to
M_{0},\>\pi_{2}(x_{1},x^{\prime})=x^{\prime}=(x_{2},\ldots,x_{n})$, are
smooth. From here, our aim is to show that $\pi_{2}$ is a local diffeomorphism
in $\partial M_{+,\varepsilon}$. To accomplish this, we note that by
definition, the vector field $\partial x_{1}$ is transversal to $\partial M$
on $\partial M_{+,\varepsilon}$. Thus, if $z_{2},\ldots,z_{n}$ are some local
coordinates in $\partial M_{+,\varepsilon}$, the functions $x_{1},z_{2},\ldots
z_{n}$ form local coordinates in $\mathbb{R}\times M_{0}$ near $z_{0}$.
Moreover, the map $x\mapsto(x_{1},x^{\prime})$ is a diffeomorphim. Thus, the
$(n-1)\times(n-1)$ matrix $\frac{\partial x_{\alpha}}{\partial z_{\beta}}$ for
$\alpha,\beta=2,\ldots,n$, which is also the differential of $\pi_{2}$, is
invertible. By the inverse function theorem, the map $\pi_{2}$ is a local
diffeomorphism in $\partial M_{+,\varepsilon}$.

Let $x\in\partial M_{+,\varepsilon}$ be an arbitrary point, and let
$\mathcal{U}\subset\partial M_{+,\varepsilon}$ be a neighborhood of $x$ such
that $\pi_{2}|_{\mathcal{U}}$ is a diffeomorphism. Then it follows from the
change of variables formula that the pullback of the surface elements satisfy
$(\pi_{2})^{\ast}(dS_{g})=J_{\pi_{2}}dV_{g_{0}}$, where $J_{\pi_{2}}$ is the
Jacobian of $\pi^{-1}_{2}$. Therefore, we have

$\int_{\mathcal{U}}|v_{s}|^{2}dS_{g}=\int_{\pi_{2}(\mathcal{U})}|v_{s}\circ\pi_{2}^{-1}|^{2}J_{\pi_{2}}dV_{g_{0}}.$

Furthermore, after possibly choosing a smaller set $\mathcal{U}$, we see that
the Jacobian $J_{\pi_{2}}$ is bounded on $\pi_{2}(\mathcal{U})\subset M_{0}$.
Therefore, we deduce from (4.35) and (4.45) that

$\int_{\mathcal{U}}|v_{s}|^{2}dS_{g}=\int_{\pi_{2}(\mathcal{U})}|v_{s}\circ\pi_{2}^{-1}|^{2}J_{\pi_{2}}dV_{g_{0}}=\mathcal{O}(1).$

Since $x\in\partial M_{+,\varepsilon}$ was arbitrarily chosen, we can choose
larger $\varepsilon$ and obtain a finite cover for $\partial
M_{+,\varepsilon}$ consisting of the sets $\mathcal{U}$ as above by shrinking
$\partial M_{+,\varepsilon}$. This leads to estimate (5.26).

Whence, estimates (5.25) and (5.26) yield

$\|v_{s}\|_{L^{2}(\Sigma_{+,\varepsilon})}+\|r_{1}\|_{L^{2}(\Sigma)}=o(h^{-1/2}),\quad
h\to 0.$ (5.27)

On the other hand, we have

$\displaystyle\bigg{(}\int_{\Sigma_{+,\varepsilon}}|\partial_{\nu}ue^{-s(\beta
t+x_{1})}|^{2}dS_{g}dt\bigg{)}^{1/2}$
$\displaystyle=\frac{1}{\sqrt{\varepsilon}}\bigg{(}\int_{\Sigma_{+,\varepsilon}}\varepsilon|\partial_{\nu}ue^{-s(\beta
t+x_{1})}|^{2}dS_{g}dt\bigg{)}^{1/2}$
$\displaystyle\leq\frac{1}{\sqrt{\varepsilon}}\bigg{(}\int_{\Sigma_{+,\varepsilon}}\partial_{\nu}\varphi|\partial_{\nu}ue^{-s(\beta
t+x_{1})}|^{2}dS_{g}dt\bigg{)}^{1/2}$
$\displaystyle\leq\frac{1}{\sqrt{\varepsilon}}\bigg{(}\int_{\Sigma_{+}}\partial_{\nu}\varphi|\partial_{\nu}ue^{-s(\beta
t+x_{1})}|^{2}dS_{g}dt\bigg{)}^{1/2},$

where we used $\Sigma_{+}=(0,T)\times\partial M_{+}^{\text{int}}$.

Using the boundary Carleman estimate (3.23) and equation (5.1), we get

$\bigg{(}\int_{\Sigma_{+}}\partial_{\nu}\varphi|\partial_{\nu}ue^{-s(\beta
t+x_{1})}|^{2}dS_{g}dt\bigg{)}^{1/2}\leq\mathcal{O}(h^{1/2})\frac{1}{\sqrt{\varepsilon}}\|e^{-s(\beta
t+x_{1})}\mathcal{L}_{c,g,a_{1},q_{1}}u\|_{L^{2}(Q)}.$

Therefore, estimates (5.24) and (5.27) yield (5.10). This completes the proof
of Lemma 5.1. ∎

###  5.2. Uniqueness of the potential

In this subsection we assume $a_{1}=a_{2}$ and prove that
$\mathcal{C}_{g,a_{1},q_{1}}=\mathcal{C}_{g,a_{2},q_{2}}$ implies
$q_{1}=q_{2}$. Our starting point is again the integral identity (5.6). When
$a_{1}-a_{2}=a=0$, this reads

$\int_{Q}qu_{2}\overline{u_{1}}dV_{g}dt=\int_{M}c^{-1}\partial_{t}u(T,x)\overline{u_{1}(T,x)}dV_{g}-\int_{\Sigma\setminus
V}\partial_{\nu}u\overline{u_{1}}dS_{g}dt.$ (5.28)

Since $a=0$, Lemma 5.1 implies that both terms on the right-hand side of
(5.28) vanish in the limit $h\to 0$. Therefore, we have

$\int_{Q}qu_{2}\overline{u_{1}}dV_{g}dt\to 0,\quad h\to 0.$

On the other hand, by substituting the CGO solutions (4.77) and (4.78) into
the left-hand side of (5.28), we get

$\int_{Q}qu_{2}\overline{u_{1}}dV_{g}dt=\int_{Q}qe^{2i\lambda(\beta
t+x_{1})}c^{-\frac{n-2}{2}}(\overline{v_{s}}w_{s}+\overline{v_{s}}r_{2}+w_{s}\overline{r_{1}}+\overline{r_{1}}r_{2})dV_{g}dt.$

It follows from estimates (4.4), (4.5), and (5.7) that

$\int_{Q}qe^{2i\lambda(\beta
t+x_{1})}c^{-\frac{n-2}{2}}(\overline{v_{s}}r_{2}+w_{s}\overline{r_{1}}+\overline{r_{1}}r_{2})dV_{g}dt=o(1),\quad
h\to 0.$

Therefore, we obtain

$\int_{Q}qe^{2i\lambda(\beta
t+x_{1})}c^{-\frac{n-2}{2}}\overline{v_{s}}w_{s}dV_{g}dt\to 0,\quad h\to 0.$

By repeating the arguments leading from (5.12) to (5.13), with the assumptions
that $q_{1},q_{2}\in C(\overline{Q})$ and $q_{1}=q_{2}$ on $\partial Q$, we
get

$\displaystyle\int_{0}^{\frac{L}{\sqrt{1-\beta^{2}}}}\int_{\mathbb{R}}\int_{\mathbb{R}}e^{2i\lambda((\beta^{2}-1)\widetilde{t}+p)-2(1-\beta^{2})\lambda
r}(cq)(\beta\widetilde{t},p-\widetilde{t},\gamma(\sqrt{1-\beta^{2}}r))$
$\displaystyle\quad\quad\quad\quad\quad\quad\times
e^{\overline{\Phi_{1}(\widetilde{t},p,r)}+\Phi_{2}(\widetilde{t},p,r)}\eta(\widetilde{t},p,r)dpd\widetilde{t}dr=0.$

Then we follow the same arguments from (5.13) onward in the proof for the
uniqueness of the damping coefficient to obtain $q_{1}=q_{2}$. This completes
the proof of Theorem 1.2.

## References

  * [1] R. A. Adams and J. J. Fournier.  Sobolev spaces.  Elsevier, 2003.
  * [2] L. Ahlfors.  Complex Analysis: An Introduction to the Theory of Analytic Functions of One Complex Variable.  McGraw-Hill Education, 3rd edition, 1979.
  * [3] S. Alexakis, A. Feizmohammadi, and L. Oksanen.  Lorentzian Calderón problem under curvature bounds.  Inventiones Mathematicae, 229:87–138, 2022.
  * [4] S. Alexakis, A. Feizmohammadi, and L. Oksanen.  Lorentzian Calderón problem near the Minkowski geometry.  Journal of European Mathematical Society, to appear, 2023.
  * [5] S. Alinhac.  Non-unicité du problème de Cauchy.  Annals of Mathematics, 117(2):77–108, 1983.
  * [6] S. Alinhac and M. Baouendi.  A non uniqueness result for operators of principal type.  Mathematische Zeitschrift, 220:561–568, 1995.
  * [7] M. Anderson, A. Katsuda, Y. Kurylev, M. Lassas, and M. Taylor.  Boundary regularity for the Ricci equation, geometric convergence, and Gel’fand’s inverse boundary problem.  Inventiones Mathematicae, 158:261–321, 2004.
  * [8] Y. Anikonov, J. Cheng, and M. Yamamoto.  A uniqueness result in an inverse hyperbolic problem with analyticity.  European Journal of Applied Mathematics, 15(5):533–543, 2004.
  * [9] V. Babich and V. Buldyrev.  Short-wavelength diffraction theory. Asymototpic methods., volume 4 of Springer Series on Wave Phenomena.  Springer (Berlin), 1991.
  * [10] V. Babich and V. Ulin.  The complex space-time ray method and quasi-photons.  Mathematical problems in the theory of wave propagation, 12, 117:5–12, 1981.
  * [11] M. Belishev.  An approach to multidimensional inverse problems for the wave equation (in Russian).  Doklady Akademii Nauk SSSR, 297(3):524–527, 1987.
  * [12] M. Belishev.  Recent progress in the boundary control method.  Inverse Problems, 23(5):R1–R67, 2007.
  * [13] M. Belishev and A. Katchalov.  Boundary control and quasiphotons in a Riemannian manifold reconstruction problem via dynamical data (in Russian).  Zap. Nauchn. Sem. Leningrad. Otdel. Mat. Inst. Steklov., 203:21–50, 1992.
  * [14] M. Belishev and Y. Kurylev.  To the reconstruction of a Riemannian manifold via its spectral data (BC-method).  Communications in Partial Differential Equations, 17(5-6):767–804, 1992.
  * [15] M. Bellassoued and I. Ben Aïcha.  An inverse problem of finding two time-dependent coefficients in second order hyperbolic equations from Dirichlet to Neumann map.  Journal of Mathematical Analysis and Applications, 475(2):1658–1684, 2019.
  * [16] M. Bellassoued and D. Dos Santos Ferreira.  Stability estimates for the anisotrpic wave equation from the Dirichlet-to-Neumann map.  Inverse Problems and Imaging, 5(4):745–773, 2011.
  * [17] M. Bellassoued, D. Jellali, and M. Yamamoto.  Lipschitz stability for a hyperbolic inverse problem by finite local boundary data.  Applicable Analysis, 85(10):1219–1243, 2006.
  * [18] M. Bellassoued and I. Rassas.  Stability estimate in the determination of a time-dependent coefficient for hyperbolic equation by partial Dirichlet-to-Neumann map.  Applicable Analysis, 98(15):2751–2782, 2019.
  * [19] M. Bellassoued and M. Yamamoto.  Determination of a coefficient in the wave equation with a single measurement.  Applicable Analysis, 87(8):901–920, 2008.
  * [20] H. Brézis.  Functional Analysis, Sobolev Spaces and Partial Differential Equations, volume 2 of Universitext.  Springer, 2011.
  * [21] M. Cekić.  The Calderón problem for connections.  Communications in Partial Differential Equations, 42(11):1781–1836, 2017.
  * [22] M. de Hoop and J. Ilmavirta.  Abel transforms with low regularity with applications to x-ray tomography on spherically symmetric manifolds.  Inverse Problems, 33:124033, 2017.
  * [23] D. Dos Santos Ferreira, C. Kenig, M. Salo, and G. Uhlmann.  Limiting Carleman weights and anisotropic inverse problems.  Inventiones Mathematicae, 178:119–171, 2009.
  * [24] D. Dos Santos Ferreira, C. Kenig, J. Sjöstrand, and G. Uhlmann.  Determining a magnetic Schrödinger operator from partial Cauchy data.  Communications in Mathematical Physics, 271(2):467–488, 2007.
  * [25] D. Dos Santos Ferreira, Y. Kurylev, M. Lassas, and M. Salo.  The Calderón problem in transversally anisotropic geometries.  Journal of European Mathematical Society, 18(11):2579–2626, 2016.
  * [26] G. Eskin.  A new approach to hyperbolic inverse problems.  Inverse Problems, 22(3):815–831, 2006.
  * [27] G. Eskin.  Inverse hyperbolic problems with time-dependent coefficients.  Communications in Partial Differential Equations, 32(11):1737–1758, 2007.
  * [28] A. Feizmohammadi, J. Ilmavirta, Y. Kian, and L. Oksanen.  Recovery of time dependent coefficients from boundary data for hyperbolic equations.  Journal of Spectral Theory, 11(3):1107–1143, 2021.
  * [29] G. Friedlander and M. Joshi.  Introduction to the Theory of Distributions.  Cambridge University Press, 1998.
  * [30] T. Helin, M. Lassas, L. Oksanen, and T. Saksala.  Correlation based passive imaging with a white noise source.  Journal de Mathématiques Pures et Appliquées, 116:132–160, 2018.
  * [31] S. O. Hussein, D. Lesnic, and M. Yamamoto.  Reconstruction of space-dependent potential and/or damping coefficients in the wave equation.  Computers & Mathematics with Applications, 6:1435–1454, 2017.
  * [32] V. Isakov.  An inverse hyperbolic problems with many boundary measurements.  Communications in Partial Differential Equations, 16:1183–1195, 1991.
  * [33] A. Katchalov and Y. Kurylev.  Multidimensional inverse problem with incomplete boundary spectral data.  Communications in Partial Differential Equations, 23(1-2):27–59, 1998.
  * [34] A. Katchalov, Y. Kurylev, and M. Lassas.  Inverse Boundary Spectral Problems, volume 123 of Monographs and Surveys in Pure and Applied Math.  Chapman & Hall/CRC, 2001.
  * [35] O. Kavian.  Lectures on parameter identification.  Three courses on partial differential equations, 4:125–162, 2003.
  * [36] C. Kenig and M. Salo.  The Calderón problem with partial data on manifolds and applications.  Analysis & PDE, 6(8):2003–2048, 2013.
  * [37] Y. Kian.  Stability of the determination of a coefficient for wave equations in an infinite waveguide.  Inverse Problems and Imaging, 8(3):713–732, 2014.
  * [38] Y. Kian.  Recovery of time-dependent damping coefficients and potentials appearing in wave equations from partial data.  SIAM Journal on Mathematical Analysis, 48(6):4021–4046, 2016.
  * [39] Y. Kian.  Stability in the determination of a time-dependent coefficient for wave equations from partial data.  Journal of Mathematical Analysis and Applications, 436(1):408–428, 2016.
  * [40] Y. Kian.  Unique determination of a time-dependent coefficient from partial data.  Annales de l’lHP (C) Analyse non linéare, 34:973–990, 2017.
  * [41] Y. Kian and L. Oksanen.  Recovery of time-dependent coefficient on Riemannian manifold for hyperbolic equations.  International Mathematics Research Notices, 2019(16):5087–5126, 2019.
  * [42] Y. Kian, L. Oksanen, and M. Morancey.  Application of the boundary coltrol method to partial data Borg-Levinson inverse spectral problem.  Mathematical Control and Related Fields, 9(2):289–312, 2019.
  * [43] K. Knudsen and M. Salo.  Determining nonsmooth first order terms from partial boundary measurements.  Inverse Problems and Imaging, 1(2):349–369, 2007.
  * [44] V. Krishnan and M. Vashisth.  An inverse problem for the relativistic Schrödinger equation with partial boundary data.  Applicable Analysis, 99(11):1889–1909, 2018.
  * [45] K. Krupchyk and G. Uhlmann.  Uniqueness on an inverse boundary problem for a magnetic Schrödinger operator with a bounded magnetic potential.  Communications in Mathematical Physics, 327(3):993–1009, 2014.
  * [46] K. Krupchyk and G. Uhlmann.  Inverse problems for magnetic Schrödinger operators in transversally anisotropic geometries.  Communications in Mathematical Physics, 361:525–582, 2018.
  * [47] I. Lasiecka, J. L. Lions, and R. Triggiani.  Nonhomogeneous boundary value problems for second order hyperbolic operators.  Journal de Mathématiques Pures et Appliquées, 65(2):149–192, 1986.
  * [48] M. Lassas and L. Oksanen.  Inverse problem for the Riemannian wave equation with Dirichlet data and Neumann data on disjoint sets.  Duke Mathematical Journal, 163(6):1071–1103, 2014.
  * [49] J. Lee.  Introduction to Smooth Manifolds, volume 218 of Graduate Texts in Mathematics.  Springer, 2nd edition, 2012.
  * [50] J. Lee.  Introduction to Riemannian Manifolds, volume 176 of Graduate Texts in Mathematics.  Springer, 2nd edition, 2018.
  * [51] R. K. Mishra and M. Vashisth.  Determining the time-dependent matrix potential in a wave equation from partial boundary data.  Applicable Analysis, 100(16):3492–3508, 2021.
  * [52] R. G. Mukhometov.  The reconstruction problem of a two-dimensional Riemannian metric, and integral geometry.  Dokl. Akad. Nauk SSSR, 232:32–35, 1977 (in Russian).
  * [53] G. Paternain, M. Salo, G. Uhlmann, and H. Zhou.  The geodesic X-ray transform with matrix weights.  American Journal of Mathematics, 141(6):1707–1750, 2019.
  * [54] Rakesh and A. G. Ramm.  Property C and an inverse problem for hyperbolic equation.  Journal of Mathematical Analysis and Applications, 156:209–219, 1991.
  * [55] Rakesh and W. Symes.  Uniqueness for an inverse problem for the wave equations.  Communications in Partial Differential Equations, 13(1):87–96, 1988.
  * [56] J. Ralston.  Approximate eigenfunctions of the Laplacian.  Journal of Differential Geometry, 12(1):87–100, 1977.
  * [57] J. Ralston.  Gaussian beams and the propagation of singularity, studies in partial differential equations.  MAA Studies in Mathematics, 23:206–248, 1982.
  * [58] A. G. Ramm and J. Sjöstrand.  An inverse problem of the wave equation.  Mathematische Zeitschrift, 206:119–130, 1991.
  * [59] R. Salazar.  Determination of time-dependent coefficients for a hyperbolic inverse problem.  Inverse Problems, 29(9):095015, 2013.
  * [60] M. Salo.  Inverse problems for nonsmooth first order perturbations of the Laplacian.  Annales Academiae Scientiarum Fennicae Mathematica Dissertationes, 139, 2004.
  * [61] S. Senapati.  Stability estimates for the relativistic Schrödinger equation from partial boundary data.  Inverse Problems, 37(1):015001, 2021.
  * [62] V. Sharafutdinov.  Integral Geometry of Tensor Fields.  Inverse and Ill-Posed Problems Series. VSP, Utrecht, 1994.
  * [63] J. Sjöstrand.  Eigenvalue distribution for non-self-adjoint operators with small multiplicative random perturbations.  Annales de la Faculté des sciences de Toulouse: Mathématiques, Serie 6, 18(4):739–795, 2009.
  * [64] P. Stefanov.  Uniqueness of the multi-dimensional inverse scattering problem for the time dependent potentials.  Mathematische Zeitschrift, 201(4):541–559, 1989.
  * [65] P. Stefanov and G. Uhlmann.  Stable determination of the hyperbolic Dirichlet-to-Neumann map for generic simple metrics.  International Mathematics Research Notices, 2005(17):1047–1061, 2005.
  * [66] D. Tataru.  Unique continuation for solutions to pde; between Hörmander’s theorem and Holmgren’s theorem.  Communications in Partial Differential Equations, 20(5-6):855–884, 1995.
  * [67] H. Whitney.  Geometric Integration Theory.  Dover Books on Mathematics. Dover Publications, 1957.
  * [68] L. Yan.  Inverse boundary problems for biharmonic operators in transversally anisotropic geometries.  SIAM Journal on Mathematical Analysis, 53(6):6617–6653, 2021.
  * [69] M. Zworski.  Semiclassical Analysis, volume 138 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, 2012.