Interpretable Asset Selection in Robust Portfolio Optimization: A Correlation Market Graph Sparsification Framework

1. Introduction

An investment portfolio is a collection of $\mathit{N}$ distinct assets, represented analytically as a vector $\mathit{w}\in {\mathbb{R}}^{\mathit{N}\times 1}$ whose entries correspond to the percentage of total wealth allocated to each asset. To avoid the risks associated with leverage and short-selling, we constrain them in $[0,1]$ and further demand that they sum up to one (i.e. ${\mathbf{1}}^T\mathit{w}=1$), satisfying the budget constraint. For every asset $\mathit{i}\in \{0,1,\dots ,\mathit{N}-1\}$, there exists a discrete-time stochastic process ${\tilde{\mathit{p}}}_i\left[n\right]$ that generates closing price data at the end of day $\mathit{n}$. Hence, we can define the daily return process of the ${\mathit{i}}^{\mathrm{th}}$ asset as

$$\widetilde{{\mathit{r}}_{\mathit{i}}}\left[\mathit{n}\right]={\displaystyle \frac{\widetilde{{\mathit{p}}_{\mathit{i}}}\left[\mathit{n}\right]-\widetilde{{\mathit{p}}_{\mathit{i}}}[\mathit{n}-1]}{\widetilde{{\mathit{p}}_{\mathit{i}}}\left[\mathit{n}\right]}}$$

(1)

Returns data has been widely adopted in the field and is generally preferred over raw price data as it typically exhibits stronger Gaussianity and serves as a metric for asset and portfolio performance [1].

A general portfolio optimization problem can be formulated as

$$\begin{array}{cc}\underset{\mathit{w}}{\min }& \left\{f\left(\mathit{w};\tilde{\mathit{r}}\right)\right\}\\ s.t.& \mathit{w}\in W\end{array}$$

(2)

where the function $f:{\mathbb{R}}^N\to \mathbb{R}$ maps a portfolio $\mathit{w}$ to a real number representing the investment risk, and is parameterized by the N-variate returns process $\tilde{\mathit{r}}\left[n\right]={\left[\begin{array}{c}{\tilde{\mathit{r}}}_0\left[n\right],\dots ,{\tilde{\mathit{r}}}_{N-1}\left[n\right]\end{array}\right]}^T$ over all portfolio assets. Note that we have also compressed the weight constraints using set notation implying that, if $\mathit{w}\in W$, then the allocation is agnostic to short selling and loans, and the budget constraint is met.

1.1. Minimum Variance Portfolio

When it comes to risk, there is no universally accepted method for quantifying it, with the choice of model often driven by complexity, robustness, or even interpretability. In 1952, H. Markowitz [2] pioneered in measuring investment risk using the mean $\mu ={\mathbb{E}}_{\tilde{\mathit{r}}}\left\{\tilde{\mathit{r}}\right\}={\left[\begin{array}{c}{\mathbb{E}}_{\tilde{\mathit{r}}}\left\{{\mathit{r}}_0\right\},{\mathbb{E}}_{\tilde{\mathit{r}}}\left\{{\mathit{r}}_1\right\},\dots ,{\mathbb{E}}_{\tilde{\mathit{r}}}\left\{{\mathit{r}}_{N-1}\right\}\end{array}\right]}^T\in {\mathbb{R}}^{N\times 1}$ and covariance $\Sigma =\mathrm{cov}\left\{\tilde{\mathit{r}}\right\}\in {\mathbb{R}}^{N\times N}$ of asset returns, which became the de-facto standard in diversification practices.

By further imposing a constraint on eq. 2 such that the portfolio expected return ${\mathit{r}}_{\mathit{p}}={\mathit{w}}^T\mu$ meets or exceeds a desirable target $\rho$ we obtain the minimum variance (MV) optimization problem as shown below

$$\begin{array}{ll}\underset{\mathit{w}}{\min }& \left\{\frac{1}{2}{\mathit{w}}^T\Sigma \mathit{w}\right\}\\ s.t.& \mathit{w}\in W\\ & {\mathbb{E}}_{\tilde{\mathit{r}}}\left({\mathit{w}}^T,\tilde{\mathit{r}}\right)\ge \rho \end{array}$$

(3)

where the cost function $f(\mathit{w};\tilde{\mathit{r}})$ is quadratic with respect to the decision variable $\mathit{w}$. Under this framework, the investor is responsible for choosing appropriate estimators for both the expected returns $\mu$ and returns covariance $\Sigma$.

The Markowitz portfolio model performs best when return distributions are not significantly skewed. However, its practical reliability is compromised due to its sensitivity to estimates of the first and second-order moments [3]. Recognizing this limitation, Markowitz himself suggested incorporating factor models to estimate the covariance matrix, particularly for large portfolios where historical data may be insufficient or unreliable [4]. These shortcomings naturally led investors to seek alternative risk metrics.

1.2. Mean Absolute Deviation Portfolio

In 1991, H. Konno and H. Yamazaki [5] proposed the use of the mean absolute deviation of portfolio returns as an alternative. They tested their work on historical data from stocks within the Japanese NIKKEI225 and TOPIX indices demonstrating that under their framework the optimization problem is more robust. This was mainly attributed to the absence of the squaring effect due to the quadratic cost in the MV problem. The MAD optimization problem was formulated as

$$\begin{array}{ll}\underset{\mathit{w}}{\min }& \left\{{\mathbb{E}}_{\tilde{\mathit{r}}}\left\{∣{\mathit{w}}^T\tilde{\mathit{r}}-{\mathit{w}}^T{\mathbb{E}}_{\tilde{\mathit{r}}}\left\{\tilde{\mathit{r}}\right\}∣\right\}\right\}\\ s.t.& \mathit{w}\in W\\ & {\mathbb{E}}_{\tilde{\mathit{r}}}\left({\mathit{w}}^T,\tilde{\mathit{r}}\right)\ge \rho \end{array}$$

(4)

where the authors arbitrarily imposed that the expected asset returns ${\mathbb{E}}_{\tilde{\mathit{r}}}\left\{\tilde{\mathit{r}}\right\}$ be estimated via the historical sample average, suggesting that

$${\mathbb{E}}_{\tilde{\mathit{r}}}\left\{\tilde{\mathit{r}}\right\}\approx {\displaystyle \frac{1}{T}}\sum _{n=0}^{T-1}\mathit{r}\left[n\right]$$

(5)

Note that we deliberately dropped the tilde symbol to indicate realisations of returns within a specified timeframe of $\mathit{T}$ samples. They also defined a term

$$\mathit{a}\left[n\right]=\left(\mathit{r}\left[n\right]-{\displaystyle \frac{1}{T}}\sum _{n=0}^{T-1}\mathit{r}\left[n\right]\right)$$

(6)

to compactly re-phrase eq. 4 as

$$\begin{array}{ll}\underset{\mathit{w}}{\min }& \left\{\frac{1}{T}{\displaystyle \sum _{n=0}^{T-1}}∣\mathit{w}T\mathit{a}\left[n\right]∣\right\}\\ s.t.& \mathit{w}\in W\\ & {\mathit{w}}^T\left(\frac{1}{T},\sum _{n=0}^{T-1},\mathit{r},\left[n\right]\right)\ge \rho \end{array}$$

(7)

The authors proceed to tractably formulate the optimisation problem in eq. 4 by splitting the absolute value and introducing $\mathit{T}$ auxiliary variables $\gamma \left[i\right]\in \mathbb{R}$ for every sample in the timeframe. The tractable formulation for the MAD optimisation problem thus becomes

$$\begin{array}{cc}\underset{\mathit{w}}{min}& \left\{{\displaystyle \frac{1}{T}}(\gamma \left[0\right]+\gamma \left[1\right]+\cdots +\gamma [T-1])\right\}\hfill \\ \left\{\gamma \right[0],\dots ,\gamma [T-1\left]\right\}& \\ \mathrm{s}.\mathrm{t}.& \mathit{w}\in W\hfill \\ & {\mathit{w}}^T\left({\displaystyle \frac{1}{T}}{\displaystyle \sum _{n=0}^{T-1}}\mathit{r}\left[n\right]\right)\ge \rho \hfill \\ & \gamma \left[n\right]+{\mathit{w}}^T\mathit{\alpha }\left[n\right]\ge 0\forall n\hfill \\ & \gamma \left[n\right]-{\mathit{w}}^T\mathit{\alpha }\left[n\right]\ge 0\forall n\hfill \end{array}$$

(8)

Konno and Yamazaki [5] also emphasize that, in contrast to the Markowitz portfolio, this involves solving a Linear Program (LP) with $2+T+T=2T+2$ constraints regardless of investment universe size.

Comparative analyses have shown that the MAD model is more efficient and robust than the MV model. G. Kasenbacher et al. [6] demonstrated that the MAD model solves portfolio optimization problems involving over 1,000 stocks much faster than the MV model by avoiding the Quadratic Programming (QP) and covariance matrix inversions required. Additionally, Y. Simaan [7] argues that the MAD model offers a more conservative measure of risk by prioritizing the spread of portfolio returns from their expected value rather than focusing solely on return variance.

While there is extensive literature on alternative risk measures, such as Keating and Shadwick’s Omega function [8] which considers higher-order moments of the returns distribution. Our primary focus, however, will be on evaluating the performance of the MAD risk metric, particularly its distributionally robust counterpart called DRMAD as proposed by D. Chen et al. [9].

2. Distributional Robustness

In real-world optimization problems, particularly in portfolio optimization, a significant challenge arises from the uncertainty in estimating parameters, as data often fails to capture the joint density governing its stochasticity, even at the first and second-order moment level [10]. This uncertainty necessitates a robust optimization framework that avoids strict assumptions about the underlying parameter distribution. Distributionally robust optimization (DRO) offers such a solution by adopting a min-max approach, where the optimization problem is framed as a game between two players [11].

In this game, one player (the minimiser) seeks to optimize the portfolio under the worst-case scenario, while the other player (the maximizer) finds this scenario by searching for the probability distribution—within an ambiguity set $\mathcal{B}$—that maximizes the expected cost. Unlike traditional robust optimization (RO), which focuses on worst-case parameter realizations within a predetermined uncertainty set U [12], DRO achieves robustness by considering the worst-case distribution itself, based on partial distributional information such as moments, support, or even distributional proximity [13]. DRO improves model conservatism and robustness against parameter uncertainty by immunizing against the worst possible data-generating distribution, rather than just the worst parameter realization, making it a valuable tool in uncertain environments.

Distributionally Robust Mean Absolute Deviation Portfolio

In 2023, Chen et al. [9] introduced a novel portfolio optimization method called DRMAD by combining the MAD model with DRO which was shown to yield superior performance in fluctuating markets as compared to classical models. In our study, their work sets the foundation with which we explore the robustness and performance of our graph-theoretic sparsification framework. According to the work of D. Chen et al. [9] the DRMAD portfolio optimization problem is defined as

$$\begin{array}{cc}\underset{\mathit{w}}{min}\hfill & \underset{\mathbb{P}}{sup}\left\{{\mathbb{E}}_{\tilde{\mathbf{r}}\sim \mathbb{P}}\left\{\left|{\mathit{w}}^T\tilde{\mathbf{r}}-{\mathit{w}}^T{\mathbb{E}}_{\tilde{\mathbf{r}}\sim \mathbb{P}}\left\{\tilde{\mathbf{r}}\right\}\right|\right\}\right\}\hfill \\ \mathrm{s}.\mathrm{t}.\hfill & {\mathbb{E}}_{\tilde{\mathbf{r}}\sim \mathbb{P}}\left\{{\mathit{w}}^T\tilde{\mathbf{r}}\right\}\ge \rho \hfill \\ \hfill & \mathbb{P}\in {\mathcal{B}}_{ϵ}\left({\widehat{\mathbb{P}}}_T\right)\hfill \\ \hfill & \mathit{w}\in W\hfill \end{array}$$

(9)

In contrast to the regular MAD problem defined in eq. 4, the distributionally robust mean absolute deviation (DRMAD) model incorporates a level of conservatism by framing the optimization task as a two-player game in which the ambiguity set, denoted as ${\mathcal{B}}_{ϵ}\left({\widehat{\mathbb{P}}}_{\mathit{T}}\right)$, is a ball that contains all distributions that are in close proximity to a distribution estimate.

To contextualise this ball-like ambiguity set, consider a finite set of realisations $\left\{\mathit{r}\right[0],\mathit{r}[1],\cdots ,\mathit{r}[\mathit{T}-1\left]\right\}$ of the returns vector $\tilde{\mathit{r}}\sim \mathbb{P}$ over a timeframe of T samples. We can estimate the unknown returns distribution ${\mathbb{P}}_{\mathit{T}}$ using a non-parametric empirical estimate defined as

$${\widehat{\mathbb{P}}}_T={\displaystyle \frac{1}{T}}\sum _{n=0}^{T-1}{\delta }_{\mathit{r}\left[n\right]}$$

(10)

where ${\delta }_{\mathit{r}\left[n\right]}$ is a Dirac measure concentrated at the point $\mathit{r}\left[n\right]$. It is reasonable to assume that the true distribution $\mathbb{P}$ is "close" to its empirical estimate ${\widehat{\mathbb{P}}}_{\mathit{T}}$ [9]. To formalise this notion of proximity, we introduce distance measures in the distribution space such as the Wasserstein distance [14].

Specificlly, given two probability distributions $\mathbb{P}$ and $\mathbb{Q}$ on the same space $\mathsf{\Xi }$, with $M\left(\mathsf{\Xi }\right)$ representing the set of distributions supported on $\mathsf{\Xi }$, their Wasserstein distance of order 1, denoted as $W_1$, is defined as

$$W_1\left(\mathbb{P},\mathbb{Q}\right)=\underset{\mathsf{\Psi }}{inf}\left\{{\int }_{\mathsf{\Xi }\times \mathsf{\Xi }}\mathrm{dist}({\tilde{\mathit{r}}}_{\mathbb{P}},{\tilde{\mathit{r}}}_{\mathbb{Q}})d\left(\mathsf{\Psi }({\tilde{\mathit{r}}}_{\mathbb{P}},{\tilde{\mathit{r}}}_{\mathbb{Q}})\right)\right\}$$

(11)

where the coupling $\mathsf{\Psi }$ is a joint distribution on $\mathsf{\Xi }\times \mathsf{\Xi }$ whose marginals are $\mathbb{P}$ and $\mathbb{Q}$, and $dist({\tilde{\mathit{r}}}_{\mathbb{P}},{\tilde{\mathit{r}}}_{\mathbb{Q}})$ is an arbitrary distance metric on the space of outcomes. Following the approach in [9], we adopt the ${\ell }_2$ norm as our distance function, i.e. $\mathrm{dist}({\tilde{\mathit{r}}}_{\mathbb{P}},{\tilde{\mathit{r}}}_{\mathbb{Q}})=\left|\right|{\tilde{\mathit{r}}}_{\mathbb{P}}-{\tilde{\mathit{r}}}_{\mathbb{Q}}{\left|\right|}^2$. This allows us to define a Wasserstein ball ${\mathbb{B}}_{ϵ}$ in the distribution space which engulfs all distributions that are within a distance $ϵ$ away from the empirical distribution ${\widehat{\mathbb{P}}}_{\mathit{T}}$. As a result the Wasserstein ambiguity set within the DRMAD formulation can be analytically defined as

$${\mathcal{B}}_{ϵ}\left({\widehat{\mathbb{P}}}_{\mathit{T}}\right)=\left\{\mathbb{P}:W_1\left(\mathbb{P},{\widehat{\mathbb{P}}}_T\right)\le ϵ\right\}.$$

(12)

As suggested by D. Kuhn and P. M. Esfahani [15,16], it is important to choose the radius hyperparameter $ϵ$ carefully to balance robustness and computational feasibility, ensuring the solution remains tractable. We also note that the radius of this ball can be determined by other metrics, such as the Kullback-Leibler divergence [17].

The DRMAD portfolio optimization problem in eq. 9 is mathematically intractable because it is formulated as an infinite-dimensional non-convex optimization problem [18]. This arises from the Wasserstein ambiguity set, where the inner maximization is performed over a distribution in an infinite-dimensional space. To overcome this challenge, D. Chen et al. [9] propose a tractable reformulation by fixing the expected portfolio value as a decision variable, which transforms the problem into a finite-dimensional convex optimization problem [9]. This reformulation results in two Linear Programs (LPs), with the optimal solution corresponding to the portfolio weights derived from the program that yields the minimum cost. According to their work, these two LPs are

$$\begin{array}{cc}\underset{\mathit{w}}{min}& \left\{t\right\}\hfill \\ t\in \mathbb{R}& \\ y_1\left[n\right],y_2\left[n\right]\in \mathbb{R}& \\ \mathrm{s}.\mathrm{t}.& t\ge {\displaystyle \frac{1}{T}}{\displaystyle \sum _{n=0}^{T-1}}{y}_1\left[n\right]+ϵ,\hfill \\ & t\ge {\displaystyle \frac{1}{T}}{\displaystyle \sum _{n=0}^{T-1}}{y}_2\left[n\right]+ϵ,\hfill \\ & y_1\left[n\right]{\mathit{w}}^T\mathit{\mu }-{\mathit{w}}^T\mathit{r}\left[n\right]-ϵ\forall n,\hfill \\ & y_1\left[n\right]\ge -{\mathit{w}}^T\mathit{\mu }+{\mathit{w}}^T\mathit{r}\left[n\right]+ϵ\forall n,\hfill \\ & y_2\left[n\right]\ge {\mathit{w}}^T\widehat{\mathit{\mu }}-{\mathit{w}}^T\mathit{r}\left[n\right]-ϵ\forall n,\hfill \\ & y_2\left[n\right]\ge -{\mathit{w}}^T\widehat{\mathit{\mu }}+{\mathit{w}}^T\mathit{r}\left[n\right]+ϵ\forall n,\hfill \\ & {\mathit{w}}^T\widehat{\mathit{\mu }}-ϵ\ge \rho ,\hfill \\ & \mathit{w}\in W\hfill \end{array}$$

(13)

$$\begin{array}{cc}\underset{\mathit{w}}{min}& \left\{t\right\}\hfill \\ t\in \mathbb{R}& \\ y_1\left[n\right],y_2\left[n\right]\in \mathbb{R}& \\ \mathrm{s}.\mathrm{t}.& t\ge {\displaystyle \frac{1}{T}}{\displaystyle \sum _{n=0}^{T-1}}{y}_1\left[n\right]+ϵ,\hfill \\ & t\ge {\displaystyle \frac{1}{T}}{\displaystyle \sum _{n=0}^{T-1}}{y}_2\left[n\right]+ϵ,\hfill \\ & y_1\left[n\right]\ge {\mathit{w}}^T\mathit{r}\left[n\right]-\rho \forall n,\hfill \\ & y_1\left[n\right]\ge -{\mathit{w}}^T\mathit{r}\left[n\right]+\rho \forall n,\hfill \\ & y_2\left[n\right]\ge {\mathit{w}}^T\widehat{\mathit{\mu }}-{\mathit{w}}^T\mathit{r}\left[n\right]-ϵ\forall n,\hfill \\ & y_2\left[n\right]\ge -{\mathit{w}}^T\widehat{\mathit{\mu }}+{\mathit{w}}^T\mathit{r}\left[n\right]+ϵ\forall n,\hfill \\ & {\mathit{w}}^T\widehat{\mathit{\mu }}-ϵ\le \rho ,\hfill \\ & {\mathit{w}}^T\widehat{\mathit{\mu }}+ϵ\ge \rho ,\hfill \\ & \mathit{w}\in W\hfill \end{array}$$

(14)

where $\mathit{w}\in {\mathbb{R}}^{N\times 1}$ represents the portfolio weights, $\widehat{\mu }={\displaystyle \frac{1}{T}}{\sum }_{n=0}^{T-1}\mathit{r}\left[n\right]$ is the sample average of the observed returns vector over a historical time window of T samples, and $y_1\left[n\right],y_2\left[n\right]\in \mathbb{R}$ are auxiliary variables for each sample. Additionally, $\rho \in \mathbb{R}$ once again represents the investor’s target returns, while the ball radius $ϵ\in [0,+\infty )$ can be interpreted as the risk appetite.

In the context of portfolio optimization, the investment universe dimension N emerges as a critical parameter. The inclusion of assets that contribute minimally or redundantly to the portfolio would complicate the optimization process by increasing the dimension of the Wasserstein ball. This would burden the search for the worst-case distribution, and would introduce unnecessary noise into the system. Additionally, a high portfolio dimensionality results in individual portfolio weights that are often too small to be of practical significance. Therefore, even though the DRMAD framework remains tractable under a large number of assets [9], it is advantageous to methodically reduce the problem dimensionality without sacrificing robustness.

3. Market Graphs

To address the inherent challenges associated with excessive portfolio dimensionality, we propose a graph-theoretic framework that sparsifies the structure of correlation market graphs, in which nodes represent assets and the edges connecting them capture historical return data correlations. This approach systematically eliminates redundant assets by reducing graph density, enabling the optimization of DRMAD portfolios within a reduced investment universe. Furthermore, this framework provides interpretability in the asset selection process prior to optimization via market graph visuals. Despite the demonstrated potential of distributionally robust optimization [9,19] and graph theory [20,21,22] in financial applications, their combined application remains unexplored in the literature. This study endeavors to fill this gap by introducing the integration of market graphs and sparsification into the DRMAD framework.

Traditionally, graphs are denoted as bold capitalised calligraphic characters, such as $\mathcal{V}$ and $\mathcal{E}$. A set of sets $\mathcal{G}=(\mathcal{V},\mathcal{E})$ is also a set and follows the same notation. Working with set notation to represent graphs, however, becomes increasingly cumbersome as the number of nodes and connections grow. To manage this, linear algebra has been established as a ubiquitous mathematical framework for modeling and analyzing graph structures. Hence, vectors are still represented by bold lowercase letters (e.g., $\mathit{r}\in {\mathbb{R}}^{N\times 1}$), and matrices are denoted by bold uppercase letters (e.g., $\mathit{A}\in {\mathbb{R}}^{N\times N}$).

In any graph, the interaction strengths between vertices u and v can be represented by scalar quantities $a_{uv}\in \mathbb{R}$. These edge weights can then be summarized in a so-called combination matrix $\mathit{A}=\left[a_{uv}\right]\in {\mathbb{R}}^{N\times N}$, where N denotes the number of vertices, i.e. the cardinality of $\mathcal{V}$ [23]. The spatial relationships among the vertices are naturally captured in the matrix structure by the arrangement of these weights. The combination matrix $\mathit{A}$ is formally constructed as

$$\mathit{A}=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1N}\\ a_{21}& a_{22}& \cdots & a_{2N}\\ ⋮& ⋮& \ddots & ⋮\\ a_{N1}& a_{N2}& \cdots & a_{NN}\end{array}\right]$$

(15)

where $a_{uv}$ indicates the interaction strength across vertices u and v.

This concise graph representation framework also allows for a logical inversion, where any square matrix can be interpreted as a graph structure, especially when its elements carry physical meaning. For example, consider a collection of N asset return realizations ${\mathit{r}}_k\left[n\right],k\in \{0,1,\cdots ,N-1\}$. The matrix $\mathit{W}$, as shown in eq. 17 below, can be defined to contain the absolute correlation coefficients $|{\rho }_{ij}|$ between all return pairs $({\mathit{r}}_i\left[n\right],{\mathit{r}}_j\left[n\right])$ in which ${\sigma }_{ij}=\mathrm{cov}({\mathit{r}}_i\left[n\right],{\mathit{r}}_j\left[n\right])$ is the sample covariance.

$$|{\rho }_{ij}|={\displaystyle \frac{|{\sigma }_{ij}|}{\sqrt{{\sigma }_{ii}{\sigma }_{jj}}}}$$

(16)

Subsequently, $\mathit{W}$ takes the form

$$\mathit{W}=\left[\begin{array}{cccc}1& |{\rho }_{12}|& \cdots & |{\rho }_{1N}|\\ |{\rho }_{21}|& 1& \cdots & |{\rho }_{2N}|\\ ⋮& ⋮& \ddots & ⋮\\ |{\rho }_{N1}|& |{\rho }_{N2}|& \cdots & 1\end{array}\right]$$

(17)

where we note that $0\le |{\rho }_{ij}|\le 1$ and $|{\rho }_{ij}|=|{\rho }_{ji}|$ by construction. To construct the correlation market graph, nodes representing assets are connected only if they are statistically correlated, with edges weighted according to the strength of this correlation. Figure 1 illustrates this using 9 randomly chosen stocks from the S&P500 index, where nodes are labeled by stock ticker symbols and edges are color-coded to qualitatively represent return data correlations.

Apart from aesthetic appeal, the market graph visual offers a geometrical interpretation of our data. For example, the strong connections between Microsoft (’MSFT’), Cisco (’CSCO’), and Intel (’INTC’) could reveal overlapping market factors such as industry, reinforcing the need for careful consideration of the investment universe prior to optimization.

Correlation graphs are typically complete, resulting in a quadratic increase in the number of edges as the number of nodes grow [23]. The maximum number of edges in an N-vertex complete graph is given by $\left({\displaystyle \genfrac{}{}{0pt}{}{N}{2}}\right)={\displaystyle \frac{1}{2}}N(N-1)$, and consequently, the combination matrix $\mathit{A}$ has no null entries, leading to market graph visuals that can become cluttered and difficult to interpret.

4. Sparsification

Graph sparsification is a process that approximates dense graphs by reducing the number of edges and nodes based on application-specific algorithms [24]. Several sparsification algorithms exist in the literature, each employing different methodologies. Zhang et. al [25] use Ricci curvature sparsification to reduce computational burden in graph neural networks (GNNs), whereas Ahmadinejad et. al [26], for example, employ a singular value approximation approach as an alternative to standard [27] (directed) graph approximations. To promote interpretability of the sparsification process, however, we employ a threshold-based sparsification algorithm on the correlation market graph $\mathit{W}$ as shown below.

Algorithm 1 Sparsification of Market Graph

Require:  Corr. matrix $\mathit{W}\in {\mathbb{R}}^{N\times N}$, corr. threshold $\tau \in [0,1]$ Ensure:  A reduced investment universe $\mathcal{S}$ of assets (set) 1: $\mathcal{S}\leftarrow \varnothing$ 2: for each $i=0$ to $N-1$ do 3:     for each $j=i$ to $N-1$ do 4:         if ${\mathit{W}}_{ij}\le \tau$ then 5:            $\mathcal{S}\leftarrow \mathcal{S}\cup \left\{{\mathrm{asset}}_i\right\}\cup \left\{{\mathrm{asset}}_j\right\}$ 6:         end if 7:     end for 8: end for 9: return$\mathcal{S}$

Given some correlation matrix $\mathit{W}\in {\mathbb{R}}^{N\times N}$ and a correlation threshold $0\le \tau \le 1$, the algorithm iterates over the upper triangular portion of $\mathit{W}$. There is no need to iterate over every pair as the correlation matrix is symmetric by construction, i.e. ${\mathit{W}}^T=\mathit{W}$. For each pair of assets, if the absolute correlation is less than or equal to the threshold $\tau$, those assets are added to the set $\mathcal{S}$ representing the new investment universe in which $\left|\mathit{S}\right|\le N$.

A visual demonstration of this sparsification algorithm is shown in Figure 2, where it has been applied on the market graph of the previous section for consistency. Moreover, we employ three correlation thresholds of 0.6, 0.4, and 0.2 to demonstrate the stock filtering capabilities of the algorithm. Note that, initially, the total number of stocks were 9, but after each sparsification procedure, only case (C) shows a reduced investment universe. In the two remaining cases the graphs were indeed sparsified but only at an edge level.

5. Numerical Example

To evaluate the performance of the DRMAD portfolio under the proposed graph sparsification scheme we used historical daily stock closing price data from the S&P500’s 100 most liquid stocks based on average trading volume. The dataset spans the period 31 DEC 2017 to 02 JAN 2024 but we split it into training and testing sets. The training data (31 DEC 2017 to 02 JAN 2019) was used to generate portfolio weights, whereas the testing dataset (02 JAN 2019 to 02 JAN 2024) was used to evaluate portfolio performance based mainly on the Sharpe ratio (SR), which in this study is defined as the ratio over the average portfolio returns to their standard deviation.

We implemented three DRMAD optimisers using the CVXPY package [28,29] with Wasserstein radii of 0, 0.1, and 1, respectively. Note that the portfolio weights from the 0 radius case are equivalent to those produced by a classic MAD portfolio as defined in eq. 8 [9]. Before deploying the optimisers, however, we also performed our graph sparsification method on the 100-asset correlation matrix using the threshold values of 0.001, 0.01, and 1. Trivially, when the correlation threshold was set to 1, the investment universe remained intact, consisting of all 100 stocks, but when using a threshold of 0.01 and 0.001 the number of assets under consideration were reduced to 32 and 9, respectively. It is common in the literature to also deploy the minimum variance (MV) portfolio (see eq. 3) whose investment universe will not be subject to sparsification. For consistency, all target returns were set to the asset-wide average training set return. Finally, we tracked the computational time needed to generate the portfolio weights under all DRMAD configurations in order to investigate the computational efficiency gain due to sparsification.

Subsequently, we applied the obtained portfolio weights to the test-set returns data spanning 5 years. To evaluate performance, we compared the growth of a $1 initial investment in each portfolio on top of the growth of the S&P500 index during the same period. Figure 3a–c display the accumulation of wealth, where each figure displays a set of DRMAD portfolios with common radius across different sparsification schemes. We also zoomed in the first six months of 2020 to observe the performance of our portfolios during unforeseen events like the breakout of the COVID-19 pandemic.

We immediately observe that all DRMAD portfolios perform the best when the Wasserstein radius is sufficiently large. In this example, a ball radius of 1 yields the best results. This can be attributed to the regularizing effect of the radius hyperparameter aiding in generalization by virtue of distributional robustness. In all of the cases however, irrespective of radius, we notice that the smaller sized DRMAD portfolios consistently outperform their full 100-stock counterparts. Indeed, there exists some correlation threshold which reduces the investment dimension without sacrificing generalization, in fact, improving it. The sparsification process should not be too intense since an investment universe that is too small could prove insufficient and potentially produce sub-optimal results.

We quantitatively summarize our findings in term of test-set return statistics and optimizer runtime on Table 1. Notice that with only 32 stocks (i.e. a third of our initial investment universe) we were able to obtain portfolios that offer higher Sharpe ratios (SR) than those with equal Wasserstein radius. The best portfolio obtained used said 32 stocks with a radius of 1 resulting in a Sharpe ratio (SR) of 0.062, surpassing the S&P500 index by 37.8% during the same 5-year period. Even with a radius of 0.1, the 32-stock portfolio achieves a superior Sharpe ratio compared to the 100-stock portfolio, even when the latter is allowed a radius ten times larger.

Finally, the best portfolio in this example highlighted with bold font on Table 1, required 75ms to compute, whereas the 100-stock equivalent needed more than triple the time.

6. Conclusions

Our study highlights the effectiveness of reducing the dimensionality of the Wasserstein DRMAD portfolio optimization method through a threshold-based sparsification procedure. By eliminating stocks with excessive statistical correlations in historical returns data, we significantly reduce the computational burden of subsequent optimization efforts and, more importantly, filter out unnecessary signal noise, thereby improving test-set generalization. Additionally, sparsifying the correlation market graph enhances robustness and promotes a faithful representation of data’s true stochastic nature, which is crucial in identifying the worst-case distribution within the Wasserstein ball centered on the empirical estimate. This novel framework, which couples DRO and sparsification, offers both visual and mathematical clarity in the decision-making process in a way that is resourceful and interpretable.