Robust and Sparse Portfolio: Optimization Models and Algorithms

Abstract

The robust and sparse portfolio selection problem is one of the most-popular and -frequently studied problems in the optimization and financial literature. By considering the uncertainty of the parameters, the goal is to construct a sparse portfolio with low volatility and decent returns, subject to other investment constraints. In this paper, we propose a new portfolio selection model, which considers the perturbation in the asset return matrix and the parameter uncertainty in the expected asset return. We define three types of stationary points of the penalty problem: the Karush–Kuhn–Tucker point, the strong Karush–Kuhn–Tucker point, and the partial minimizer. We analyze the relationship between these stationary points and the local/global minimizer of the penalty model under mild conditions. We design a penalty alternating-direction method to obtain the solutions. Compared with several existing portfolio models on seven real-world datasets, extensive numerical experiments demonstrate the robustness and effectiveness of our model in generating lower volatility.

1. Introduction

In 1952, Harry M. Markowitz [1] published the classic “Portfolio Selection” in The Journal of Finance, which ushered in a new era of financial mathematical analysis. Markowitz pointed out that investors who care about return and risk should hold portfolios located at the efficient boundary of mean-variance, which is the famous mean-variance portfolio (MVP) selection model. Since then, many portfolio selection strategies have been proposed by referring to the MVP and its variants. However, MVPs exhibit instability due to estimation errors in the input parameters [2], especially in large-scale conditions. The instability means that the solution obtained under sample fluctuation may be optimal for a given sample, but it is not optimal from the perspective of risk. For more comments on this model, we refer to [3,4,5,6] and the references therein.

This paper focuses attention on sample fluctuations and parameter uncertainty in the portfolio selection problem. We now review some relevant methods for the parameter uncertainty. Among various approaches, the attractive one is the robust portfolio (RP), which corresponds to a robust optimization, since it does not use any information about the probability distribution of the uncertain parameters. RP we considered is a conservative approach that minimizes the loss function within an uncertainty set and then solves the problem under the worst-case scenario. In the last two decades, robust portfolio selection problems have gained the increasing interest of researches. These researches constructed well-known optimal portfolios from the perspective of robust optimization [7,8,9,10]. In this way, Goldfarb and Iyengar [11] formulated and solved RP problems. They introduced the uncertainty structures for the input parameters, then they showed that the RP problems corresponding to the second-order cone programs and these uncertainty structures correspond to confidence regions employed to estimate the market parameters. Given the uncertainty in the mean and covariance matrix of the asset return, Lobo and Boyd [12] computed the maximum risk of a portfolio in a numerically efficient way. They proved that this is a semi-definite programming problem and is readily solved by interior-point methods for convex optimization. Min et al. [13] proposed the hybrid RP models under ellipsoidal uncertainty sets, and they considered both the best-case and the worst-case counterparts. Won and Kim [14] considered RP problems involving a trade-off between the worst-case utility and the worst-case regret, or the largest difference between the best utility achievable under the model and that achieved by a given portfolio. They showed that the entire optimal trade-off curve can be found via solving a series of semi-definite programs under the ellipsoidal uncertainty model. Some research works [15,16] concentrated on the application of robust optimization on basic mean-variance, mean value-at-risk (mean-VaR), and mean conditional-value-at-risk (mean-CVaR) problems, but did not consider variants of the problem like robust index tracking, robust and sparse portfolio selection problems, and so on. More relevant works can be found in [17,18,19,20] and the references therein.

RPs have a wide range of applications, among these, one essential step is the construction of uncertainty sets. Two types of uncertainty sets are widely used, namely the box uncertainty set and the ellipsoidal uncertainty set. Tütüncü and Koenig [21] used symmetric box uncertainty sets defined as $U_{\mu }=\{\mu \in {\mathbb{R}}^n|{\mu }^L\le \mu \le {\mu }^U\}$ and $U_{\mathrm{\Sigma }}=\{\mathrm{\Sigma }\in {\mathbb{R}}^{n\times n}|{\mathrm{\Sigma }}^L\le \mathrm{\Sigma }\le {\mathrm{\Sigma }}^U,\mathrm{\Sigma }⪰0\}$, where ${\mu }^L\in {\mathbb{R}}^n$ and ${\mu }^U\in {\mathbb{R}}^n$ are the lower and upper bounds of mean vector $\mu$, ${\mathrm{\Sigma }}^L\in {\mathbb{R}}^{n\times n}$ and ${\mathrm{\Sigma }}^U\in {\mathbb{R}}^{n\times n}$ are the lower and the upper bounds of the covariance matrix $\mathrm{\Sigma }$, respectively, and $\mathrm{\Sigma }$ is positive semi-definite. Khodamoradi et al. [22] used box uncertainty sets for a cardinal-constrained mean-variance portfolio problem which allows short selling. Swain and Ojha [10] analyzed the robust version of the mean-variance portfolio problem and mean-semi-variance portfolio problem with box uncertainty sets. Alternatively, Fabozzi et al. [23] defined an ellipsoidal uncertainty set for the expected asset return as $U_{\mu }=\left\{\mu \right|(\mu -\overline{\mu }){\mathrm{\Sigma }}^{-1}{(\mu -\overline{\mu })}^{\top }\le {ϵ}^2\}$, where $\overline{\mu }$ is the nominal asset return and ${ϵ}^2$ is a small scalar, which controls the size of the uncertainty set. However, they did not consider the uncertainty of the covariance matrix, thus the solution was robust only against perturbations in the asset return vector. Pıinar [24] developed a multi-period robust mean-variance portfolio problem with an ellipsoidal uncertainty set while allowing short selling. As we all know, the estimation error is more sensitive to the mean vector than the covariance matrix. On the other hand, dealing with the uncertainty in the covariance matrix is more complicated than dealing with the uncertainty set of the mean vector. Thus, in this paper, we consider two types of uncertainty sets for the mean vector.

Financial data have some remarkable features, such as multicollinearity and a heavy tail. Therefore, the perturbations of these data should not be underestimated. By referring to Brodie et al. [2], who transferred the MVP into a Lasso-type portfolio, we consider the perturbations in the asset return matrix and design its uncertainty set. In addition, from the perspective of transaction costs and administrative expenses, more assets are not always better. Therefore, it is also necessary to consider sparsity when constructing a portfolio [25,26,27]. After these discussions, a natural question follows: How do we find better RPs that not only reduce the undesired impact of parameter uncertainty, but also improve sparsity and reduce cost?

Following the above considerations, this paper proposes a sparsity constrained robust portfolio optimization model with parameter uncertainty and data perturbation. Specifically, we consider the perturbation in the asset return matrix and the parameter uncertainty in the expected asset return. By using the equivalence of robustness and regularization, the Lasso-type objective function can be converted into the sum of a square root and the ${\ell }_1$ norm. We consider two kinds of uncertainty sets: the box uncertainty set and the ellipsoidal uncertainty set. For its penalty model, we define three types of stationary points: the Karush–Kuhn–Tucker (KKT) point, the strong KKT point, and the partial minimizer. Under mild constraint qualification (CQ), we prove that any local minimizer of the penalty model is a KKT point. Moreover, the global minimizer of the penalty model is proven to be a partial minimizer and, then, a stronger KKT point under Slater’s CQ. Finally, a penalty alternating direction method is proposed to obtain a portfolio, and its convergence is established. We confirm the effectiveness of our approach by comparing with nine widely studied portfolio models on seven real-world data sets. The numerical results show that the portfolios we proposed have less volatility, that is less risk. Moreover, our portfolio strategies can yield higher Sharpe ratios when the appropriate parameters are selected.

This paper is organized as follows. Some notations and preliminaries used in this paper are given in the next section. The model of robust and sparse portfolios and the analysis of their optimization theory are stated in Section 3. Two types of uncertainty sets of mean vectors are presented in Section 4. The optimization algorithm named the penalty alternating direction method is established in Section 5. Extensive numerical experiments are conducted in Section 6. Conclusions are drawn in Section 7.

2. Notations and Preliminary

We use $\mathbb{R}$ and ${\mathbb{R}}^n$ and ${\mathbb{R}}^{m\times n}$ to denote the set of real numbers and the n-dimensional and $m\times n$-dimensional Euclidean space. We use boldfaced small letters to denote vectors, e.g., $\mathbf{w}\in {\mathbb{R}}^n$ is a column vector with n elements $w_i$, $i=1,\dots ,n$. The transpose of $\mathbf{w}$ is denoted as ${\mathbf{w}}^T$, which is a row vector. In particular, ${\mathbf{1}}_n$ is the vector of all ones of size n. For a vector $\mathbf{a}\in {\mathbb{R}}^n$, we define its absolute value vector by $\left|\mathbf{a}\right|:=\left(\right|a_1|,\cdots ,|a_n\left|\right)$. We use capital letters to denote matrices, e.g., $A\in {\mathbb{R}}^{m\times n}$ and $a_{ij}$ denote the $(i,j)$-th entry of A. Given an index $\Gamma \subset \{1,\dots ,n\}$, ${\mathbf{a}}_{\Gamma }$ denotes the sub-vector of $\mathbf{a}$. We write the Euclidean norm of $\mathbf{w}$ by ${\parallel \mathbf{w}\parallel }_2$, the ${\ell }_1$ norm by ${\parallel \mathbf{w}\parallel }_1$, and the infinity norm by ${\parallel \mathbf{w}\parallel }_{\infty }$. For two vectors $\mathbf{a}\in {\mathbb{R}}^n$ and $\mathbf{b}\in {\mathbb{R}}^n$, $\langle \mathbf{a},\mathbf{b}\rangle$ denotes the standard inner product.

We now provide some existing results of optimization that are crucial for the theory of this paper. For the convenience of expression, we define the following convex programming:

$$\begin{array}{cc}\hfill \underset{\mathbf{x}\in {\mathbb{R}}^n}{min}& f\left(\mathbf{x}\right),\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& g_i\left(\mathbf{x}\right)\ge 0,i=1,\dots ,m,\mathbf{x}\in \mathsf{\Omega },\hfill \end{array}$$

(1)

where $\mathsf{\Omega }$ is a nonempty convex set, f is a convex function, and the $g_i\left(\mathbf{x}\right)$s are concave functions. For problem (1), Slater’s CQ builds a bridge between its solution and the KKT point (the point satisfying the conditions in Theorem 1).

Definition 1

([28], Definition 4.17). Slater’s CQ holds in problem (1) if there exists $\mathbf{u}\in \mathsf{\Omega }$ such that $g_i\left(u\right)>0$ for all $i=1,\dots ,m$.

Theorem 1

([28], Theorem 4.18). Suppose that Slater’s CQ holds in problem (1). Then, ${\mathbf{x}}^{*}$ is an optimal solution to problem (1) if and only if there exist non-negative Lagrange multipliers $({\lambda }_1,\dots ,{\lambda }_m)\in {\mathbb{R}}^m$ such that

$$0\in \partial f\left({\mathbf{x}}^{*}\right)-\sum _{i=1}^m{\lambda }_i\partial g_i\left({\mathbf{x}}^{*}\right)+N({\mathbf{x}}^{*};\mathsf{\Omega })$$

and ${\lambda }_i\partial g_i\left({\mathbf{x}}^{*}\right)=0$ for all $i=1,\dots ,m$, where $\partial f\left({\mathbf{x}}^{*}\right)$ denotes the classical sub-differential set ([28], Definition 2.30) of f at ${\mathbf{x}}^{*}$ and $N({\mathbf{x}}^{*};\#)$ denotes the classical normal cone ([28], Definition 2.9) of # at ${\mathbf{x}}^{*}$.

We also introduce some crucial terminologies and results for sparsity nonlinear programming:

$$\begin{array}{cc}\hfill \underset{\mathbf{x}\in {\mathbb{R}}^n}{min}& f\left(\mathbf{x}\right)\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& g_i\left(\mathbf{x}\right)\ge 0,i=1,\dots ,m,\hfill \\ \hfill & h_j\left(\mathbf{x}\right)=0,j=1,\dots ,l,\hfill \\ \hfill & {\parallel \mathbf{x}\parallel }_0\le s,\hfill \end{array}$$

(2)

where f is a convex function and g and h are continuously differentiable. A restricted linear independence constraint qualification (R-LICQ) used for sparsity nonlinear programming (2) was defined by [29] as follows.

Definition 2

([29], Definition 2.4). We say that the R-LICQ holds at ${\mathbf{x}}^{*}$, where ${\mathbf{x}}^{*}$ is feasible for the problem (2):

• When ${\parallel {\mathbf{x}}^{*}\parallel }_0=s$, $\nabla g_i\left({\mathbf{x}}^{*}\right)$, $i\in I\left({\mathbf{x}}^{*}\right)$, $\nabla h_j\left({\mathbf{x}}^{*}\right)$, $j=1,\dots ,l,$ are linearly independent.
• When ${\parallel {\mathbf{x}}^{*}\parallel }_0<s$, ${\nabla }_{{\mathrm{\Gamma }}^{*}}{g}_i\left({\mathbf{x}}^{*}\right)$, $i\in I\left({\mathbf{x}}^{*}\right)$, ${\nabla }_{{\Gamma }^{*}}{h}_j\left({\mathbf{x}}^{*}\right)$, $j=1,\dots ,l,$ are linearly independent.

Based on the R-LICQ, the following decomposition result holds.

Theorem 2

([29], Proposition 2.5). Let ${\mathbf{x}}^{*}$ be a feasible point of problem (2) and the R-LICQ hold at ${\mathbf{x}}^{*}$. Then,

$$\widehat{N}({\mathbf{x}}^{*};S\cap Q)=\widehat{N}({\mathbf{x}}^{*};S)+\widehat{N}({\mathbf{x}}^{*};Q),$$

where $S:=\{\mathbf{x}:\parallel \mathbf{x}{\parallel }_0\le s\}$, $Q:=\{\mathbf{x}:g_i\left(\mathbf{x}\right)\ge 0,i=1,\dots ,m,h_j\left(\mathbf{x}\right)=0,j=1,\dots ,l\}$, and $\widehat{N}({\mathbf{x}}^{*};\#)$ denotes the Frechét normal cone ([30], Definition 6.3) of # at ${\mathbf{x}}^{*}$, which degenerates into the classical norm cone described in Theorem 1 if # is a convex set.

For the partial problem (10) of the portfolio model (6) in Section 3.1, the R-LICQ holds automatically at ${\mathbf{x}}^{*}$, where $m=0$, $l=1$, and $h\left(\mathbf{x}\right):={\mathbf{1}}^T\mathbf{x}-1$. Next, we establish the relationship between the local minimizer of problem (2) and its KKT point (the point satisfying the KKT system in Theorem 3).

Theorem 3.

Suppose that ${\mathbf{x}}^{*}$ is a local minimizer of problem (2) and the R-LICQ holds at ${\mathbf{x}}^{*}$. Then, there exist non-negative Lagrange multipliers $({\lambda }_1^{*},\dots ,{\lambda }_m^{*})\in {\mathbb{R}}_{+}^m$ and $({\mu }_1^{*},\dots ,{\mu }_m^{*})\in {\mathbb{R}}^l$ such that

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill & 0\in \partial f\left({\mathbf{x}}^{*}\right)-\sum _{i=1}^m{\lambda }_i\partial g_i\left({\mathbf{x}}^{*}\right)+\sum _{j=1}^l{\lambda }_i\partial h_j\left({\mathbf{x}}^{*}\right)+\widehat{N}({\mathbf{x}}^{*};S),\hfill \\ \hfill & g_i\left(\mathbf{x}\right)\ge 0,{\lambda }_i{g}_i\left(\mathbf{x}\right)=0,i=1,\dots ,m,\hfill \\ \hfill & h_j\left(\mathbf{x}\right)=0,j=1,\dots ,l,\hfill \\ \hfill & {\parallel \mathbf{x}\parallel }_0\le s.\hfill \end{array}\right.\end{array}$$

(3)

Proof. 

It follows from Theorem 6.12 of [30] that

$$0\in \partial f\left({\mathbf{x}}^{*}\right)+\widehat{N}({\mathbf{x}}^{*};S\cap Q).$$

Combing Theorem 2 with the proof of Theorem 3.2 of [29], this result holds.    □

This result is different from Theorem [29]. We allow the objective function of problem (2) to be non-differentiable. The analysis process of this result is completely consistent with that of Theorem [29].

3. Model and Optimization Theory

In this section, we first propose a robust and sparse portfolio model (4) with an uncertainty set constraint and a sparsity constraint. For the convenience of the numerical calculation, we consider its ${\ell }_1$ norm penalization variant (6). We define three types of stationary points of the penalization variant: the KKT point, the strong KKT point, and the partial minimizer. The relationships of these stationary points and the local/global minimizer of the penalization problem (6) are established in Section 3.2.

3.1. Robust and Sparse Portfolio Model

Consider n risky assets, denoting the asset return at period t by ${\mathbf{r}}_t={(r_1,...,r_n)}^{\top }\in {\mathbb{R}}^n$. The expected return vector of different assets is denoted by $E\left({\mathbf{r}}_t\right)=\mu$, and the covariance matrix is denoted by $E\left[({\mathbf{r}}_t-\mu ){({\mathbf{r}}_t-\mu )}^{\top }\right]=V$. In the traditional Markowitz portfolio selection problem, the portfolio construction is based on the trade-off between risk and return. For a given level of acceptable portfolio return $\rho ={\mathbf{w}}^{\top }\mu$, the mean-variance optimization can be formulated as

$$\begin{array}{c}\hfill \underset{\mathbf{w}\in {\mathbb{R}}^n}{min}\frac{1}{2}{\mathbf{w}}^{\top }V\mathbf{w},s.t.{\mathbf{w}}^{\top }\mu =\rho ,{\mathbf{w}}^{\top }{\mathbf{1}}_n=1,\end{array}$$

and its aim is to find a portfolio that has minimal risk for a given expected return. A significant model that has been developed from the Markowitz model is the Lasso-type portfolio proposed by Brodie et al. [2], which is given as:

$$\begin{array}{c}\hfill \underset{\mathbf{w}\in {\mathbb{R}}^n}{min}\frac{1}{T}\parallel \rho {\mathbf{1}}_T{-R\mathbf{w}\parallel }_2^2+\alpha {\parallel \mathbf{w}\parallel }_1,s.t.{\mathbf{w}}^{\top }\overline{\mu }=\rho ,{\mathbf{w}}^{\top }{\mathbf{1}}_n=1,\end{array}$$

where $\overline{\mu }=\frac{1}{T}{\sum }_{t=1}^T{\mathbf{r}}_t$, $\alpha$ is the penalty parameter, and $R\in {\mathbb{R}}^{T\times n}$ is the asset return matrix. Brodie et al. [2] confirmed that the ${\ell }_1$ norm can produce a sparse portfolio, and this method can stabilize the problem. In this paper, we start with the square root Lasso-type portfolio, while adding more consideration about the perturbation in asset return matrix R and the parameter uncertainty in $\mu$. We propose the following robust and sparse portfolio selection model:

$$\begin{array}{cc}\hfill \underset{\mathbf{w}\in {\mathbb{R}}^n}{min}\underset{\Delta \in U_0}{max}& \parallel \rho {\mathbf{1}}_T{-(R+\Delta )\mathbf{w}\parallel }_2\hfill \\ \hfill s.t.& \underset{\mu \in U}{min}{\mathbf{w}}^T\mu \ge \rho ,\hfill \\ \hfill & {\mathbf{w}}^T{\mathbf{1}}_n=1,\mathbf{w}\ge 0,\hfill \\ \hfill & {\parallel \mathbf{w}\parallel }_0\le s,\hfill \end{array}$$

(4)

where $\Delta$ is the data perturbation matrix and $U_0=\{\Delta \in {\mathbb{R}}^{T\times n}:\parallel {\Delta }_i{\parallel }_2\le \alpha ,\forall i\in \{1,...,n\}\}$. The uncertainty set of the asset return is denoted by U, and we will discuss two selections of U in the last section.

In [31] (Chapter 2), they showed the equivalence of robustness and regularization. Specifically, they precisely characterized the conditions on the model of uncertainty and loss function under which robustness is equivalent to regularization for linear regression.

Definition 3.

Let $g:{\mathbb{R}}^T\to \mathbb{R}$ and $h:{\mathbb{R}}^n\to \mathbb{R}$ be the norm, then the induced norm ${\parallel ·\parallel }_{(h,g)}$ is defined as

$${\parallel \Delta \parallel }_{(h,g)}=\underset{\mathbf{w}\in {\mathbb{R}}^n}{max}\frac{g(\Delta \mathbf{w})}{h\left(\mathbf{w}\right)}.$$

Theorem 4

([31], Chapter 2). If $r,q\in [1,\infty ]$, then

$$\underset{\mathbf{w}}{min}\underset{\Delta \in U_{({\ell }_q,{\ell }_r)}}{max}{\parallel \mathbf{y}-(R+\Delta )\mathbf{w}\parallel }_r=\underset{\mathbf{w}}{min}{\parallel \mathbf{y}-R\mathbf{w}\parallel }_r+\alpha {\parallel \mathbf{w}\parallel }_q,$$

where $U_{({\ell }_q,{\ell }_r)}={\{\Delta :\parallel \Delta \parallel }_{({\ell }_q,{\ell }_r)}\le \alpha \}$. Moreover, if $U_0=\{\Delta :\parallel {\Delta }_i{\parallel }_2\le \alpha ,\forall i\in \{1,...,n\}\}$, then $U_{({\ell }_1,{\ell }_2)}=U_0$, and this implies

$$\underset{\mathbf{w}}{min}\underset{\parallel {\Delta }_i{\parallel }_2\le \alpha }{max}{\parallel \mathbf{y}-(R+\Delta )\mathbf{w}\parallel }_2=\underset{\mathbf{w}}{min}{\parallel \mathbf{y}-R\mathbf{w}\parallel }_2+\alpha {\parallel \mathbf{w}\parallel }_1.$$

From the relationship of the robustness and the regularization, problem (4) can be rewritten as

$$\begin{array}{cc}\hfill \underset{\mathbf{w}\in {\mathbb{R}}^n}{min}& \parallel \rho {\mathbf{1}}_T{-R\mathbf{w}\parallel }_2+\alpha {\parallel \mathbf{w}\parallel }_1\hfill \\ \hfill s.t.& \underset{\mu \in U}{min}{\mathbf{w}}^T\mu \ge \rho ,\hfill \\ \hfill & {\mathbf{w}}^T{\mathbf{1}}_n=1,\mathbf{w}\ge 0,\hfill \\ \hfill & {\parallel \mathbf{w}\parallel }_0\le s.\hfill \end{array}$$

(5)

Under this transformation, the problem (5) actually enjoys robustness. We plan to use an alternating penalty method to solve problem (5). To ensure the implementation of the alternating penalty method, we add a copy constraint $\mathbf{w}=\mathbf{v}$ to the problem (5) and, then, move it to the objective function by means of the ${\ell }_1$ norm penalty, then the penalization formulation is

$$\begin{array}{cc}\hfill \underset{\mathbf{w},\mathbf{v}\in {\mathbb{R}}^n}{min}& f(\mathbf{w},\mathbf{v}):=\parallel \rho {\mathbf{1}}_T{-R\mathbf{w}\parallel }_2+{\alpha \parallel \mathbf{w}\parallel }_1+\beta {\parallel \mathbf{w}-\mathbf{v}\parallel }_1\hfill \\ \hfill s.t.& \mathbf{w}\in {\mathsf{\Omega }}_1:=\left\{\mathbf{w}\right|\underset{\mu \in U}{min}{\mathbf{w}}^T\mu \ge \rho ,\mathbf{w}\ge 0\}\hfill \\ \hfill & \mathbf{v}\in {\mathsf{\Omega }}_2:=\left\{\mathbf{v}\right|{\mathbf{v}}^T{\mathbf{1}}_n={1,\parallel \mathbf{v}\parallel }_0\le s\}.\hfill \end{array}$$

(6)

We conduct its optimality analysis in the next subsection.

3.2. Optimization Theory

We now analyze the optimality of the penalization problem (6). Obviously, the objective function of the problem (6) is a lower semi-continuous and coercive function. Theorem 5 in the next subsection provides the existence of optimal solutions. The selection of the uncertainty set U is discussed in Section 4.

This subsection provides a few theoretical results of the problem (6) including the existence of the solution and three classes of the first-order necessary optimal condition.

Theorem 5.

For any given $\alpha \in {\mathbb{R}}_{+}$ and $\beta \in {\mathbb{R}}_{+}$, the optimal solutions of the problem (6) can be attained.

Proof. 

It is clear that f is a proper, closed, and coercive function and ${\mathsf{\Omega }}_1\times {\mathsf{\Omega }}_2$ is a nonempty closed set satisfying ${\mathsf{\Omega }}_1\times {\mathsf{\Omega }}_2\cap \mathrm{dom}\left(f\right)\ne \varnothing$. It follows from Theorem 2.14 of [32] that this theorem holds.    □

We now define a class of KKT points of the problem (6). For the convenience of expression and the generality of optimality, we write ${min}_{\mu \in U}{\mathbf{w}}^T\mu \ge \rho$ as $g\left(\mathbf{w}\right)\ge 0$ and suppose that g is a concave function and is not necessarily differentiable. Indeed, the quadratic uncertainty set and the absolute uncertainty set introduced in Section 4 satisfy these terminologies.

Definition 4.

The point $({\mathbf{w}}^{*},{\mathbf{v}}^{*})\in {\mathsf{\Omega }}_1\times {\mathsf{\Omega }}_2$ is called a KKT point of the problem (6), if there exist Lagrange multipliers ${\lambda }_1^{*}\in {\mathbb{R}}_{+}$ and ${\lambda }_2^{*}\in \mathbb{R}$ such that the following system holds:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill & 0\in {\partial }_{\mathbf{w}}f({\mathbf{w}}^{*},{\mathbf{v}}^{*})-{\lambda }_1^{*}{\partial }_{\mathbf{w}}g\left({\mathbf{w}}^{*}\right)+N({\mathbf{w}}^{*};{\mathbb{R}}_{+}),\hfill \\ \hfill & 0\in {\partial }_{\mathbf{v}}f({\mathbf{w}}^{*},{\mathbf{v}}^{*})-{\lambda }_2^{*}\mathbf{1}+\widehat{N}({\mathbf{v}}^{*};S),\hfill \\ \hfill & g\left({\mathbf{w}}^{*}\right)\ge 0,{\lambda }^{*}g\left({\mathbf{w}}^{*}\right)=0,\hfill \\ \hfill & {\mathbf{1}}^T{\mathbf{v}}^{*}=1,{\parallel {\mathbf{v}}^{*}\parallel }_0\le s.\hfill \end{array}\right.\end{array}$$

(7)

Although the functions corresponding to the quadratic uncertainty set and the absolute uncertainty set introduced in Section 4 are all concave and may both be non-differentiable and Slater’s CQ automatically holds for both functions, we still considered Slater’s CQ as a condition of Theorem 6 for the sake of generality. Moreover, it is stated in Section 2 that the R-LICQ of ${\mathsf{\Omega }}_2$ holds at every point. Then, only under the condition that Slater’s CQ holds, the relationship between the local minimizer of the problem (6) and the KKT point of the problem (6) can be established.

Theorem 6.

Let $({\mathbf{w}}^{*},{\mathbf{v}}^{*})\in {\mathsf{\Omega }}_1\times {\mathsf{\Omega }}_2$ be a local minimizer of the problem (6). If Slater’s CQ holds on ${\mathsf{\Omega }}_1$, then it is a KKT point of the problem (6).

Proof. 

On the one hand, since $({\mathbf{w}}^{*},{\mathbf{v}}^{*})\in {\mathsf{\Omega }}_1\times {\mathsf{\Omega }}_2$ is a local minimizer of the problem (6), ${\mathbf{w}}^{*}$ is a local minimizer of the following optimization:

$$\begin{array}{cc}\hfill \underset{\mathbf{w}\in {\mathbb{R}}^n}{min}& f(\mathbf{w},{\mathbf{v}}^{*})=\parallel \rho {\mathbf{1}}_T{-R\mathbf{w}\parallel }_2+{\alpha \parallel \mathbf{w}\parallel }_1+\beta {\parallel \mathbf{w}-{\mathbf{v}}^{*}\parallel }_1\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathbf{w}\in {\mathsf{\Omega }}_1.\hfill \end{array}$$

(8)

Notice that $f(\mathbf{w},{\mathbf{v}}^{*})$ is a convex function about $\mathbf{w}$ and ${\mathsf{\Omega }}_1$ is a convex set. Then, problem (8) is a convex optimization. Since Slater’s CQ holds on ${\mathsf{\Omega }}_1$, it follows from Theorem 1 that there exists a Lagrange multiplier ${\lambda }_1^{*}\in {\mathbb{R}}_{+}$ such that

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill & 0\in {\partial }_{\mathbf{w}}f({\mathbf{w}}^{*},{\mathbf{v}}^{*})+{\lambda }_1^{*}{\partial }_{\mathbf{w}}g\left({\mathbf{w}}^{*}\right)+N({\mathbf{w}}^{*};{\mathbb{R}}_{+}),\hfill \\ \hfill & g\left({\mathbf{w}}^{*}\right)\ge 0,{\lambda }_1^{*}g\left({\mathbf{w}}^{*}\right)=0.\hfill \end{array}\right.\end{array}$$

(9)

On the other hand, ${\mathbf{v}}^{*}$ is a local minimizer of the following optimization:

$$\begin{array}{cc}\hfill \underset{\mathbf{v}\in {\mathbb{R}}^n}{min}& f({\mathbf{w}}^{*},\mathbf{v})=\parallel \rho {\mathbf{1}}_T-R{\mathbf{w}}^{*}{\parallel }_2+\alpha \parallel {\mathbf{w}}^{*}{\parallel }_1+\beta {\parallel {\mathbf{w}}^{*}-\mathbf{v}\parallel }_1\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathbf{v}\in {\mathsf{\Omega }}_2.\hfill \end{array}$$

(10)

Since the R-LICQ of ${\mathsf{\Omega }}_2$ holds at every point, it follows from Theorem 3 that there exists a Lagrange multiplier ${\lambda }_2^{*}\in \mathbb{R}$ such that

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill & 0\in {\partial }_{\mathbf{v}}f({\mathbf{w}}^{*},{\mathbf{v}}^{*})+{\lambda }_2^{*}\mathbf{1}+\widehat{N}({\mathbf{v}}^{*};S),\hfill \\ \hfill & {\mathbf{1}}^T{\mathbf{v}}^{*}=1,{\parallel {\mathbf{v}}^{*}\parallel }_0\le s.\hfill \end{array}\right.\end{array}$$

(11)

Combing the system (9) and (11), this theorem holds.    □

Again, problem (10) can be simply written as

$$\begin{array}{cc}\hfill \underset{\mathbf{v}\in {\mathbb{R}}^n}{min}& \parallel {\mathbf{w}}^{*}{-\mathbf{v}\parallel }_1\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathbf{v}\in {\mathsf{\Omega }}_2,\hfill \end{array}$$

and it has a closed-form solution; see [33], i.e.,

$$\begin{array}{c}\hfill v_i^{*}=\left\{\begin{array}{cc}\frac{w_i^{*}}{{\left({\mathbf{w}}_s^{*}\right)}^T{\mathbf{1}}_s},\hfill & \mathrm{if}i\in I_s^{*}\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.\end{array}$$

(12)

where $I_s^{*}:=\left\{i\right|w_{\langle 1\rangle }^{*}\ge \dots \ge w_{\langle s\rangle }^{*}\}$ and $w_{\langle i\rangle }^{*}$ denotes the i-th largest absolute value among the n elements of ${\mathbf{w}}^{*}$. Thus, we can define a class of strong KKT points of the problem (6) as follows.

Definition 5.

The point $({\mathbf{w}}^{*},{\mathbf{v}}^{*})\in {\mathsf{\Omega }}_1\times {\mathsf{\Omega }}_2$ is called a strong KKT point of the problem (6), if there exists a Lagrange multiplier ${\lambda }^{*}\in {\mathbb{R}}_{+}$ such that the following system holds:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill & 0\in {\partial }_{\mathbf{w}}f({\mathbf{w}}^{*},{\mathbf{v}}^{*})+{\partial }_{\mathbf{w}}g\left({\mathbf{w}}^{*}\right)+N({\mathbf{w}}^{*};{\mathbb{R}}_{+}),\hfill \\ \hfill & g\left({\mathbf{w}}^{*}\right)\ge 0,{\lambda }^{*}g\left({\mathbf{w}}^{*}\right)=0,\hfill \\ \hfill & v_i^{*}=\left\{\begin{array}{cc}\frac{w_i^{*}}{{\left({\mathbf{w}}_s^{*}\right)}^T{\mathbf{1}}_s},\hfill & \mathrm{if}i\in I_s^{*}\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}\right.\end{array}$$

(13)

It is easy to prove that, if $({\mathbf{w}}^{*},{\mathbf{v}}^{*})$ is a strong KKT point of the problem (6), then it is a KKT point of the problem (6). The following result provides the relationship between the global minimizer of the problem (6) and the strong KKT point of the problem (6).

Theorem 7.

Let $({\mathbf{w}}^{*},{\mathbf{v}}^{*})\in {\mathsf{\Omega }}_1\times {\mathsf{\Omega }}_2$ be a global minimizer of the problem (6). If Slater’s CQ holds on ${\mathsf{\Omega }}_1$ at $w^{*}$, then it is a strong KKT point of the problem (6).

Proof. 

The part of ${\mathbf{w}}^{*}$ in (13) follows from (7). We only need to discuss the part of ${\mathbf{v}}^{*}$ in (13). Since ${\mathbf{v}}^{*}$ is the global minimizer of (10), it follows from (12) that the part of ${\mathbf{v}}^{*}$ in (13) holds.    □

Note that the local minimizer of the problem (6) cannot be guaranteed to be a strong KKT point.

Finally, we introduce the third stationary point of the problem (6), which is called the partial minimizer.

Definition 6.

The point $({\mathbf{w}}^{*},{\mathbf{v}}^{*})\in {\mathsf{\Omega }}_1\times {\mathsf{\Omega }}_2$ is called a partial minimizer of the problem (6), if it satisfies

$$f({\mathbf{w}}^{*},{\mathbf{v}}^{*})\le f(\mathbf{w},{\mathbf{v}}^{*}),\forall \mathbf{w}\in {\mathsf{\Omega }}_1,f({\mathbf{w}}^{*},{\mathbf{v}}^{*})\le f({\mathbf{w}}^{*},\mathbf{v}),\forall \mathbf{v}\in {\mathsf{\Omega }}_2.$$

Clearly, any global minimizer of the problem (6) is a partial minimizer. Moreover, on the one hand, the partial problem (8) is a convex optimization, and Slater’s CQ ensures that its KKT point and global minimizer are consistent. On the other hand, the partial problem (10) has a closed-form solution. Thus, the equivalence relationship between the KKT point of the problem (6) and the partial minimizer of the problem (6) can be established under Slater’s CQ.

Theorem 8.

Let $({\mathbf{w}}^{*},{\mathbf{v}}^{*})\in {\mathsf{\Omega }}_1\times {\mathsf{\Omega }}_2$ be a feasible point of the problem (6). Suppose that Slater’s CQ holds on ${\mathsf{\Omega }}_1$. Then, $({\mathbf{w}}^{*},{\mathbf{v}}^{*})$ is a partial minimizer of the problem (6) if and only if $({\mathbf{w}}^{*},{\mathbf{v}}^{*})$ is a strong KKT point of the problem (6).

Proof. 

Suppose that $({\mathbf{w}}^{*},{\mathbf{v}}^{*})$ is a strong KKT point of the problem (6), then

$$0\in {\partial }_{\mathbf{w}}f({\mathbf{w}}^{*},{\mathbf{v}}^{*})+{\partial }_{\mathbf{w}}g\left({\mathbf{w}}^{*}\right)+N({\mathbf{w}}^{*};{\mathbb{R}}_{+}),g\left({\mathbf{w}}^{*}\right)\ge 0,\mathrm{and}{\lambda }_1^{*}g\left({\mathbf{w}}^{*}\right)=0.$$

Since Slater’s CQ holds at ${\mathbf{w}}^{*}$, ${\mathbf{w}}^{*}$ is a global minimizer of the problem (8). Then, we have that $f({\mathbf{w}}^{*},{\mathbf{v}}^{*})\le f(\mathbf{w},{\mathbf{v}}^{*})$$\forall \mathbf{w}\in {\mathsf{\Omega }}_1$. Moreover, it follows from the definition of the strong KKT point of the problem (6) that ${\mathbf{v}}^{*}$ is a global minimizer of the problem (10). Then, we have that $f({\mathbf{w}}^{*},{\mathbf{v}}^{*})\le f({\mathbf{w}}^{*},\mathbf{v})$, $\forall \mathbf{v}\in {\mathsf{\Omega }}_2$. Thus, $({\mathbf{w}}^{*},{\mathbf{v}}^{*})$ is a partial minimum of the problem (6). The opposite conclusion clearly holds.    □

4. The Uncertainty Set U

In SubSection 3.2, we rewrite the uncertainty set constraint as $g\left(w\right)\ge 0$, where g is a generalized concave function and is not necessarily differentiable. This section introduces two mainstream formulations for the uncertainty set in asset mean return vector $\mu$ (see [34]), which corresponds to the quadratic uncertainty set and the absolute uncertainty set, respectively.

4.1. The Quadratic Uncertainty Set

The first one is the quadratic formulation, $U=\left\{\mu \right|{(\mu -\overline{\mu })}^T\mathsf{\Omega }(\mu -\overline{\mu })\le {\kappa }^2\}$, where $\overline{\mu }$ is the nominal expected return and $\kappa$ is the error. Assume that asset returns are independent and identically distributed and $\mu -\overline{\mu }$ follows a normal distribution with mean value $\mathbf{0}$ and covariance matrix $\mathsf{\Omega }$, where $\mathsf{\Omega }$ is the covariance matrix of errors in the expected asset return. In Yin et al. [35], they discussed the choice of uncertainty matrix $\mathsf{\Omega }$ in the quadratic uncertainty set and proposed the selection criteria. In the quadratic uncertainty case, ${min}_{\mu \in U}{\mathbf{w}}^T\mu$ in problem (5) is equivalent to the following problem:

$$\underset{\mu \in U}{max}{\mathbf{w}}^T\overline{\mu }-{\mathbf{w}}^T\mu .$$

Solving the above problem, we obtain

$$\mu =\overline{\mu }-\sqrt{\frac{{\kappa }^2}{{\mathbf{w}}^T\mathsf{\Omega }\mathbf{w}}}\mathsf{\Omega }\mathbf{w}.$$

Then, the problem (5) is rewritten as

$$\begin{array}{cc}\hfill \underset{\mathbf{w},\mathbf{v}\in {\mathbb{R}}^n}{min}& \parallel \rho {\mathbf{1}}_T{-R\mathbf{w}\parallel }_2+\alpha {\parallel \mathbf{w}\parallel }_1\hfill \\ \hfill s.t.& {\overline{\mu }}^T\mathbf{w}-\kappa \sqrt{{\mathbf{w}}^T\mathsf{\Omega }\mathbf{w}}\ge \rho ,\hfill \\ \hfill & {\mathbf{w}}^T{\mathbf{1}}_n=1,\mathbf{w}\ge 0,\hfill \\ \hfill & {\parallel \mathbf{w}\parallel }_0\le s.\hfill \end{array}$$

(14)

Here, $g\left(\mathbf{w}\right)={\overline{\mu }}^T\mathbf{w}-\kappa \sqrt{{\mathbf{w}}^T\mathsf{\Omega }\mathbf{w}}-\rho$. The penalization form of problem (14) can be rewritten as

$$\begin{array}{cc}\hfill \underset{\mathbf{w},\mathbf{v}\in {\mathbb{R}}^n}{min}& \parallel \rho {\mathbf{1}}_T{-R\mathbf{w}\parallel }_2+{\alpha \parallel \mathbf{w}\parallel }_1+\beta {\parallel \mathbf{w}-\mathbf{v}\parallel }_1\hfill \\ \hfill s.t.& \mathbf{w}\in {\mathsf{\Omega }}_1=\left\{\mathbf{w}\right|\kappa \parallel \sqrt{\mathsf{\Omega }}\mathbf{w}{\parallel }_2\le {\overline{\mu }}^T\mathbf{w}-\rho ,\mathbf{w}\ge 0\}\hfill \\ \hfill & \mathbf{v}\in {\mathsf{\Omega }}_2=\left\{\mathbf{v}\right|{\mathbf{v}}^T{\mathbf{1}}_n={1,\parallel \mathbf{v}\parallel }_0\le s\}.\hfill \end{array}$$

According to the proposition of Yin et al. [35], we choose $\mathsf{\Omega }\propto diag\left(V\right)$. By using this uncertainty matrix, it is expected to reduce the sensitivity to the inputs, as well as keep the original volatility unchanged.

4.2. The Absolute Uncertainty Set

Fabozzi et al. [23] used the absolute uncertainty set in mean returns that ask that the sum of absolute spreads between estimated and possible mean returns should not be too large. The absolute formulation is $U=\left\{\mu \right|{\sum }_i|{\mu }_i-{\overline{\mu }}_i|\le \frac{\kappa \overline{\sigma }}{\sqrt{T}}\}$. In this case,

$${\mu }^T\mathbf{w}-{\overline{\mu }}^T\mathbf{w}\ge -\sum _i|{\mu }_i-{\overline{\mu }}_i|max(|{\mathbf{w}}_i\left|\right)\ge -\frac{\kappa \overline{\sigma }}{\sqrt{T}}max\left(\right|{\mathbf{w}}_i\left|\right);$$

thus,

$${\mu }^T\mathbf{w}\ge {\overline{\mu }}^T\mathbf{w}-\frac{\kappa \overline{\sigma }}{\sqrt{T}}max\left(\right|{\mathbf{w}}_i\left|\right).$$

Then, the problem (5) is equivalent to

$$\begin{array}{cc}\hfill \underset{\mathbf{w},\mathbf{v}\in {\mathbb{R}}^n}{min}& \parallel \rho {\mathbf{1}}_T{-R\mathbf{w}\parallel }_2+\alpha {\parallel \mathbf{w}\parallel }_1\hfill \\ \hfill s.t.& {\overline{\mu }}^T\mathbf{w}-\frac{\kappa \overline{\sigma }}{\sqrt{T}}max\left(\right|{\mathbf{w}}_i\left|\right)\ge \rho ,\hfill \\ \hfill & {\mathbf{w}}^T{\mathbf{1}}_n=1,\mathbf{w}\ge 0,\hfill \\ \hfill & {\parallel \mathbf{w}\parallel }_0\le s.\hfill \end{array}$$

(15)

Here, $g\left(\mathbf{w}\right)={\overline{\mu }}^T\mathbf{w}-\kappa \overline{\sigma }max\left(\right|{\mathbf{w}}_i\left|\right)/\sqrt{T}-\rho$. Similarly, the penalization from of problem (15) can be written as

$$\begin{array}{cc}\hfill \underset{\mathbf{w},\mathbf{v}\in {\mathbb{R}}^n}{min}& \parallel \rho {\mathbf{1}}_T{-R\mathbf{w}\parallel }_2+{\alpha \parallel \mathbf{w}\parallel }_1+\beta {\parallel \mathbf{w}-\mathbf{v}\parallel }_1\hfill \\ \hfill s.t.& \mathbf{w}\in {\mathsf{\Omega }}_1=\left\{\mathbf{w}\right||w_i|\le \frac{\sqrt{T}}{\kappa \overline{\sigma }}({\mathbf{w}}^T\overline{\mu }-\rho ),i=1,...,n,\mathbf{w}\ge 0\}\hfill \\ \hfill & \mathbf{v}\in {\mathsf{\Omega }}_2=\left\{\mathbf{v}\right|{\mathbf{v}}^T{\mathbf{1}}_n={1,\parallel \mathbf{v}\parallel }_0\le s\}.\hfill \end{array}$$

5. Optimization

This section introduces a penalty alternating direction method (PADM) to solve problem (5).

5.1. Alternating Direction Methods

We first discuss the optimization of the problem (6). Due to the complexity of this problem, alternating direction methods (ADMs) can be used to solve this problem. The framework of ADMs is described as follows:

Next, we state the general convergence result of Algorithm 1, and one can refer to Geissler et al. [36] for a proof (Theorem 8) and for further details about this method.

Algorithm 1 ADM: Alternating Direction Method.

1: Set the problem parameters: $\alpha ,\kappa ,\rho ,T>0$, asset return matrix $R\in {\mathbb{R}}^{T\times n}$, and nominal expected return vector $\overline{\mu }\in {\mathbb{R}}^n$. Initialize $\epsilon >0$, $({\mathbf{w}}^0$, ${\mathbf{v}}^0)$ and penalty parameter $\beta >0$. Set the iteration index $k:=0,1,...$. 2: Compute $${\mathbf{w}}^{k+1}\in arg\underset{\mathbf{w}}{min}\{f(\mathbf{w},{\mathbf{v}}^k):\mathbf{w}\in {\mathsf{\Omega }}_1\},$$ (16) and $${\mathbf{v}}^{k+1}\in arg\underset{\mathbf{v}}{min}\{f({\mathbf{w}}^{k+1},\mathbf{v}):\mathbf{v}\in {\mathsf{\Omega }}_2\}.$$ (17) 3: If $\frac{\parallel {\mathbf{w}}^{k+1}-{\mathbf{w}}^k{\parallel }_2+{\parallel {\mathbf{v}}^{k+1}-{\mathbf{v}}^k\parallel }_2}{\parallel {\mathbf{w}}^k{\parallel }_2+{\parallel {\mathbf{v}}^k\parallel }_2}\le \epsilon$, then stop with $({\mathbf{w}}^k,{\mathbf{v}}^k)$ being an output point of (6).

Theorem 9.

Let $\left\{({\mathbf{w}}^k,{\mathbf{v}}^k)\right\}$ be a sequence generated by Algorithm 1. Then, the following holds:

(a)

$\left\{({\mathbf{w}}^k,{\mathbf{v}}^k)\right\}$ is bounded.

(b)

Any limiting point $\left\{({\mathbf{w}}^{*},{\mathbf{v}}^{*})\right\}$ of $\left\{({\mathbf{w}}^k,{\mathbf{v}}^k)\right\}$ is a partial minimizer of the problem (6).

(c)

If Slater’s CQ holds on ${\mathsf{\Omega }}_1$, the limiting point of $\left\{({\mathbf{w}}^k,{\mathbf{v}}^k)\right\}$ is also a strong KKT of the problem (6).

Proof. 

(a) It follows from Algorithm 1 that

$$f({\mathbf{w}}^{k+1},{\mathbf{v}}^{k+1})\le f({\mathbf{w}}^{k+1},{\mathbf{v}}^k)\le f({\mathbf{w}}^k,{\mathbf{v}}^k).$$

Since f is a coercive function, then the level set of f is bounded. Thus, $\left\{({\mathbf{w}}^k,{\mathbf{v}}^k)\right\}$ is bounded.

(b) Clearly, $\left\{f({\mathbf{w}}^k,{\mathbf{v}}^k)\right\}$ is a decreasing sequence and $f({\mathbf{w}}^k,{\mathbf{v}}^k)\ge 0$, then there exists a value $f^{*}$ such that ${lim}_{k\to \infty }f({\mathbf{w}}^k,{\mathbf{v}}^k)=f^{*}$. Suppose that $\left\{({\mathbf{w}}^{*},{\mathbf{v}}^{*})\right\}$ is a limiting point of $\left\{({\mathbf{w}}^k,{\mathbf{v}}^k)\right\}$. Then, there exists a sequence $\left\{k_j\right\}$ such that ${lim}_{j\to \infty }{k}_j=\infty$, ${lim}_{j\to \infty }({\mathbf{w}}^{k_j},{\mathbf{v}}^{k_j})=({\mathbf{w}}^{*},{\mathbf{v}}^{*})$ and ${lim}_{j\to \infty }f({\mathbf{w}}^{k_j},{\mathbf{v}}^{k_j})=f({\mathbf{w}}^{*},{\mathbf{v}}^{*})=f^{*}$, where the last equality holds since ${lim}_{k\to \infty }f({\mathbf{w}}^k,{\mathbf{v}}^k)=f^{*}$. Without loss of generality, let ${lim}_{k\to \infty }({\mathbf{w}}^k,{\mathbf{v}}^k)=({\mathbf{w}}^{*},{\mathbf{v}}^{*})$. It follows from Algorithm 1 that

$${\mathbf{w}}^{k+1}\in arg\underset{\mathbf{w}\in {\mathsf{\Omega }}_1}{min}f(\mathbf{w},{\mathbf{v}}^k).$$

Since f is continuous with respect to $\mathbf{w}$, then taking $k\to \infty$, we have that

$${\mathbf{w}}^{*}\in arg\underset{\mathbf{w}\in {\mathsf{\Omega }}_1}{min}f(\mathbf{w},{\mathbf{v}}^{*}).$$

Similarly,

$${\mathbf{v}}^{*}\in arg\underset{\mathbf{v}\in {\mathsf{\Omega }}_2}{min}f({\mathbf{w}}^{*},\mathbf{v}).$$

Thus, the limiting point $\left\{({\mathbf{w}}^{*},{\mathbf{v}}^{*})\right\}$ of $\left\{({\mathbf{w}}^k,{\mathbf{v}}^k)\right\}$ is a partial minimizer of the problem (6).

(c) Under Slater’s CQ, the partial minimizer of the problem (6) is a strong KKT point of the problem (6) and the opposite also holds. Thus, this result holds.    □

5.2. The Optimization for the Partial Problem (8)

We now discuss the optimization of the partial problem (8) at the k-th iteration of the ADM. Some non-exact penalty methods and smoothing methods can be used to solve this problem. Here, we obtain ${\mathbf{w}}^{k+1}$ by solving the following optimization:

$$\begin{array}{cc}\hfill \underset{\mathbf{w}\in {\mathbb{R}}^n}{min}& f(\mathbf{w},{\mathbf{v}}^k)=\parallel \rho {\mathbf{1}}_T{-R\mathbf{w}\parallel }_2+{\alpha \parallel \mathbf{w}\parallel }_1+\beta \parallel \mathbf{w}-{\mathbf{v}}^k{\parallel }_1+\gamma \left|g{\left(\mathbf{w}\right)}_{-}\right|\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathbf{w}\ge 0,\hfill \end{array}$$

(18)

where $\gamma >0$ is a penalty parameter. Let

$${\psi }_{\mu }\left(t\right)=\left\{\begin{array}{cc}\left|t\right|,\hfill & \left|t\right|\ge \mu ,\hfill \\ \frac{t^2}{2\mu }+\frac{\mu }{2},\hfill & \left|t\right|<\mu ,\hfill \end{array}\right.{\varphi }_{\mu }\left(t\right)=\frac{1}{2}(t+\sqrt{t^2+\mu }),$$

where $\mu >0$ is a smoothing parameter. Then, a class of the smoothing optimization of problem (18) can be given as follows:

$$\begin{array}{cc}\hfill \underset{\mathbf{w}\in {\mathbb{R}}^n}{min}& f_{\mu }(\mathbf{w},{\mathbf{v}}^k)=\sqrt{\parallel \rho {\mathbf{1}}_T{-R\mathbf{w}\parallel }_2^2+\mu }+\hfill \\ \hfill & \alpha \sum _{i=1}^n{\psi }_{\mu }\left(w_i\right)+\beta \sum _{i=1}^n{\psi }_{\mu }(w_i-v_i^k)+\gamma {\varphi }_{\mu }(-g\left(\mathbf{w}\right))\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathbf{w}\ge 0.\hfill \end{array}$$

(19)

The projection gradient method (PGM) can be used to solve the problem (19), and its iteration formula is

$${\mathbf{w}}^{j+1,k}=P_{{\mathbb{R}}_{+}}({\mathbf{w}}^{j,k}+\eta \nabla f_{\mu }({\mathbf{w}}^{j,k},{\mathbf{v}}^k)),$$

where $\eta >0$ denotes the step length at the j-th iteration of the PGM at the k-th iteration of the ADM and $P_{\#}\left(t\right)$ denotes the projection point of t onto #. The framework of the above method is called the penalty projection gradient method (PPGM) and can be described in Algorithm 2.

Algorithm 2 PPGM: Penalty Projection Gradient Method.

1: Set the problem parameters: $\alpha ,\kappa ,\rho ,T,\beta ,{\gamma }_{max}>0$, asset return matrix $R\in {\mathbb{R}}^{T\times n}$, and nominal expected return vector $\overline{\mu }\in {\mathbb{R}}^n$. Initialize penalty parameters ${\gamma }^0>0$, $\tau >1$. Set the iteration index $j=0,1,\cdots$. 2: Computing ${\mathbf{w}}^{j+1,k}=P_{{\mathbb{R}}_{+}}({\mathbf{w}}^{j-1,k}+\eta \nabla f_{\mu }(\mathbf{w},{\mathbf{v}}^k))$. 3: If ${\mathbf{w}}^{j+1,k}$ satisfies $\frac{\parallel {\mathbf{w}}^{j+1,k}-{\mathbf{w}}^{j,k}{\parallel }_2}{\parallel {\mathbf{w}}^{j,k}{\parallel }_2}\le \epsilon$ and $g\left({\mathbf{w}}^{j,k}\right)\ge -\epsilon$, then stop and ${\mathbf{w}}^k={\mathbf{w}}^{j,k}$. 4: If $\frac{\parallel {\mathbf{w}}^{j+1,k}-{\mathbf{w}}^{j,k}{\parallel }_2}{\parallel {\mathbf{w}}^{j,k}{\parallel }_2}\le \epsilon$ and $g\left(w^{j,k}\right)<-\epsilon$, then choose new penalty parameter $\gamma =min\{\tau \gamma ,{\gamma }_{max}\}$. Otherwise, return to Step 2.

5.3. Penalty Alternating Direction Method

At the end of this section, we describe the PADM for the general problem (5). At iteration l, set the value of penalty parameter ${\beta }^l$ and obtain $({\mathbf{w}}^l,{\mathbf{v}}^l)$ by the ADM with ${\beta }^l$. If the inequality $\parallel {\mathbf{w}}^l-{\mathbf{v}}^l{\parallel }_1\le tol$ holds, where $tol$ is a small positive constant, we stop with a feasible solution of problem (6). Otherwise, the penalty parameter ${\beta }^l$ is updated to ${\beta }^{l+1}$. In this way, the PADM generates a sequence of the partial minimizer of problem (6) with ${\beta }^l$. The framework of the PADM is formally stated in Algorithm 3.

Algorithm 3 PADM: Penalty Alternating Direction Method.

1: Set the problem parameters: $\alpha ,\kappa ,\rho ,T,{\beta }_{max}>0$, asset return matrix $R\in {\mathbb{R}}^{T\times n}$, and nominal expected return vector $\overline{\mu }\in {\mathbb{R}}^n$. Initialize penalty parameters ${\beta }^0>0$, $\tau >1$. Set the iteration index $l=0,1,\cdots$. 2: Obtain $({\mathbf{w}}^l,{\mathbf{v}}^l)$ by the ADM with ${\beta }^l$. 3: If $({\mathbf{w}}^l,{\mathbf{v}}^l)$ satisfies $\parallel {\mathbf{w}}^l-{\mathbf{v}}^l{\parallel }_1\le tol$, then stop with $({\mathbf{w}}^l,{\mathbf{v}}^l)$. Otherwise, choose new penalty parameter ${\beta }^{l+1}=min\{\tau {\beta }^l,{\beta }_{max}\}$, and return to Step 2.

6. Numerical Results

This section shows extensive numerical experiments. In Section 6.1, we first present six real data sets, explain some existing models to be compared with and describe the performance measures to be used. In Section 6.2, we demonstrate that our methods lead to robust and sparse portfolios. In Section 6.3, we compare nine popular portfolios in terms of out-of-sample (OOS) performance measures. Finally, in Section 6.4, we show the cumulative return of different portfolio strategies. All of our computations are conducted in the Matlab R2019a environment, on a PC with an Intel(R) Core(TM) i5-7200U CPU (2.50 GHz, 4 CPUs) and 4G RAM processors.

6.1. Models of Comparison, Data, and Performance Measures

(a) Eleven portfolio models compared. We compare the OOS performance of 11 portfolio models across six real data sets of weekly and monthly returns. Those models are well studied, and we divide them into four groups, which are summarized in

Table 1. List of portfolio strategies considered.

Group	Model	Abbr.	Refer.
(1)	Robust and sparse portfolios with
	quadratic uncertainty set	RSQ	this paper
	absolute uncertainty set	RSA	this paper
(2)	Some well-studied portfolio strategies with
	${\ell }_1$ regularization	L1	Brodie et al. [2]
	${\ell }_{1,2}$ regularization	L12	Zhao et al. [37]
	Elastic Net regularization	EN	Yen and Yen [38]
	upper and lower bound	Box	Behr et al. [39]
(3)	Benchmarks’ portfolio strategies with
	short-sales constrained	SC	Jagannathan and Ma [5]
	short-sales unconstrained	SU	Jagannathan and Ma [5]
	equally weighted (1/N) portfolio	EW	DeMiguel et al. [40]
(4)	Shrinkage of covariance
	sample covariance and identity matrix	SCID	Olivier and Wolf [41]
	sample covariance and 1-factor matrix	SC1F	Olivier and Wolf [3]

. The first group is the robust and sparse portfolio strategies developed in this paper. The second group includes some well-studied portfolio strategies. The third group includes three benchmark portfolio strategies. The last group consists of two portfolios that use the shrinkage technique to estimate the covariance matrix.

(b) Seven data sets tested.

Table 2. Information of the seven real data sets.

#	Data Sets	Stocks	Time Period	Source	Frequency
1	DJIA	29	01/10/2017–30/10/2022	Yahoo finance	Weekly
2	NASDAQ	95	01/10/2017–30/10/2022	Yahoo finance	Weekly
3	SP500	336	01/10/2017–30/10/2022	Yahoo finance	Weekly
4	Russell2000	1340	01/10/2017–30/10/2022	Yahoo finance	Weekly
5	Russell3000	2166	01/10/2017–30/10/2022	Yahoo finance	Weekly
6	SP100	71	01/10/2017–30/10/2022	Yahoo finance	Weekly
7	FF100	100	11/1999–06/2022	K.French	Monthly

lists some real-world data sets: DJIA [42], NASDAQ [43], S&P [44,45], Russell2000 [46], Russell3000 [47], and FF100 [38]. All the data are obtained from Yahoo finance (https://finance.yahoo.com/, accessed on 10 January 2023) and Ken French’s website (https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html, accessed on 10 January 2023). In all cases, we remove those assets that have missing values.

(c) Measuring the OOS performance and its setup. We largely follow the “rolling-window” procedures in [2,37] to conduct our comparison. Let T be the length of a data set and $\tau$ be the window length (e.g., $\tau =120$) used to construct the optimal portfolio by a model. In each period $(t+1)$, $t=\tau ,...,T-1$, we compute different portfolios over the previous $\tau$ periods. We then compute the OOS return in the $(t+1)$-th period based on the obtained portfolio. We repeat this procedure until we reach the end of the data set. In this way, we will obtain a series of $(T-\tau )$ portfolio vectors for each model listed in Table 1. To make it precise, let ${\mathbf{w}}_t^s$ be the optimal portfolio obtained by the portfolio strategy s over the date from $t-\tau +1,\dots ,t$. The OOS return in the $t+1$ period is computed as ${\mathbf{r}}_{t+1}^s={{\mathbf{w}}_t^s}^{⊺}{\mathbf{r}}_{t+1}$, where ${\mathbf{r}}_{t+1}$ is the return in the $(t+1)$-th period. Thus, we obtain a time series of $(T-\tau -1)$ periods OOS returns for all strategies. Note that we use the traditional “rolling-window” procedures for the numerical analysis, and some new methods could provide new ideas for the analysis of portfolio selection problems, see [48].

The OOS performance of each portfolio strategy is assessed by using four quantities: (i) the OOS portfolio variance (${\widehat{\sigma }}^2$), (ii) the OOS portfolio Sharpe ratio ($\widehat{SR}$), (iii) portfolio turnover ($TURN$), and (iv) the average short positions ($ASP$). The specific definitions can be found in DeMiguel et al. [6], Yen and Yen [38], and Zhao et al. [37]. We evaluate the cumulative return (CR). The CR of a portfolio scores the total payoffs that are yielded by the investment strategy across the investment periods without considering any risk or cost, see Shen et al. [49]. We also consider some quantities studied in [38] on the profiles of the portfolio weights: PAP represents the proportion of active positions and PZP is the proportion of zero positions, respectively, defined as $PA{P}_t=\frac{|S_t^1|}{N}$, $PZ{P}_t=\frac{|S_t^0|}{N}$, where $S_t^1=\{i:w_{i,t}\ne 0\}$ and $S_t^0=\{i:w_{i,t}=0\}$.

6.2. Robust and Sparse Portfolio

This section shows the weight of robust and sparse portfolios. We use the DJIA data set and the sparse levels $s_1=\lceil 30\%n\rceil$ and $s_2=\lceil 50\%n\rceil$. The parameter $\alpha =\beta =10\lambda$, and the value of $\lambda$ varies from ${10}^{-2}$ to ${10}^1$.

Figure 1 shows the portfolio weights, PAP, and PZP. The two plots in the top panel correspond to a robust portfolio under the quadratic uncertainty set, and the sparsity is $s_1$. The two plots in the bottom panel correspond to a robust portfolio under the absolute uncertainty set, and the sparsity is $s_2$. With the increase of penalty parameter $\lambda$, the portfolio weights tended to be sparse. The PAP and PZP indicate that we can obtain sparse portfolios that satisfy the specified sparsity.

Figure 2 shows the sparse portfolio. We use four different data sets. The sparsity level on DJIA is $s_1=\lceil 30\%n\rceil$, on NASDAQ and FF100 is $s_2=\lceil 10\%n\rceil$, and on Russell2000 is $s_3=\lceil 1\%n\rceil$. We solve the robust portfolio under the quadratic uncertainty set to show the results. We obtain the portfolio with the specified sparsity and the distribution of different asset weight values.

6.3. Out-of-Sample Performance

The Sharpe ratio considers return and risk at the same time; it is a comprehensive measurement for us to observe the performance of a portfolio. Thus, we first test the Sharpe ratio of different portfolio strategies. We use the SP100 data set. The parameter $\alpha =2\beta$ and the value of $\beta$ varies from 10 to 10${}^{1.5}$. The sparsity level $s_1=\lceil 15\%n\rceil$ and $s_2=\lceil 5\%n\rceil$.

By comparing with two benchmark portfolios, Figure 3 shows that the RSQ and RSA can produce a higher Sharpe ratio when choosing a suitable penalty parameter.

Table 3. Portfolio out-of-sample variance (${\widehat{\sigma }}^2$) (${(\%)}^2$), Sharpe ratio ($\widehat{SR}$), turnover ($TURN$), and the average short positions ($ASP$).

		DJIA	NASDAQ	SP500	Russell2000	Russell3000	FF100
		n = 29	n = 95	n = 336	n = 1340	n = 2166	n = 100
	var	5.3741	5.9065	11.5604	12.1813	12.7837	17.2215
RSA	SR	0.0703	0.1807	0.1008	0.0063	0.0267	0.3215
	TURN	0.1194	0.1721	0.1430	0.1033	0.1021	0.1698
	ASP	−1.11e-18	2.22e-18	0	−3.08e-18	−4.01e-18	−6.28e-18
	var	5.3949	5.8633	12.5936	12.7041	12.7702	17.2097
RSQ	SR	0.0737	0.1889	0.0863	0.0059	0.0258	0.3216
	TURN	0.1063	0.1711	0.1334	0.1112	0.1003	0.1698
	ASP	−6.66e-18	−3.33e-18	−1.43e-17	−9.25e-18	−9.25e-18	0
	var	9.5931	7.1442	32.1907	14.8927	13.7879	18.3918
L12	SR	0.0860	0.1768	0.1000	0.0072	0.0295	0.2971
	TURN	0.0220	0.0286	0.0369	0.0675	0.0605	0.0296
	ASP	−0.0216	−0.0209	−0.0240	−0.0057	−0.0050	−0.0276
	var	8.0373	6.7078	13.1228	14.2679	13.0166	16.4937
L1	SR	0.0864	0.1831	0.0414	0.0042	0.0279	0.3279
	TURN	0.1063	0.0677	0.0369	0.1250	0.1136	0.0629
	ASP	−0.0036	−0.0012	−0.0124	−0.0039	0.0020	−0.0195
	var	8.5871	6.5362	12.9449	14.2074	13.0817	16.5823
EN	SR	0.0911	0.1790	0.0418	0.0029	0.0402	0.3181
	TURN	0.0256	0.0472	0.0408	0.1208	0.1181	0.0510
	ASP	−0.0257	−0.0198	−0.0272	0.0015	0.0011	−0.0257
	var	9.9181	8.7749	3.59e+02	7.3636	7.0004	34.5878
BOX	SR	0.0250	0.0768	−0.1204	−0.0495	−0.0023	0.6244
	TURN	0.8322	1.7146	3.2709	0.2820	0.2973	5.1895
	ASP	1.1172	2.7850	6.5137	0.5744	0.5886	8.4327
	var	10.0132	8.0089	86.1291	24.3110	16.9741	21.2942
SC	SR	0.0871	0.1891	0.1202	0.0233	0.0399	0.2812
	TURN	0.0388	0.0313	0.0411	0.0512	0.0481	0.0247
	ASP	1.38e-16	1.23e-16	−1.52e-16	3.12e-16	−4.19e-16	0
	var	9.9181	15.0721	5.91e+03	12.1721	9.9999	17.2369
SU	SR	0.0250	0.1800	−0.1270	−0.0482	−0.0167	0.4624
	TURN	0.8322	2.7947	3.9150	0.3919	0.3722	5.2429
	ASP	1.1172	2.5580	5.8870	0.5542	0.5494	6.3456
	var	11.1339	8.1780	89.8252	21.3370	19.5742	22.2346
EW	SR	0.0701	0.1790	0.1151	0.0194	0.0337	0.2922
	TURN	0.0208	0.0253	0.0400	0.0429	0.0381	0.0252
	ASP	1.13e-16	1.13e-16	−1.12e-16	4.52e-16	−3.39e-16	0
	var	6.9600	6.9782	2.38e+03	7.3676	7.0304	16.8727
SCID	SR	0.0294	0.1097	−0.1259	−0.0497	−0.0025	0.6740
	TURN	0.4115	0.8589	0.9940	0.2803	0.2884	1.5393
	ASP	0.6362	1.6735	2.8275	0.5727	0.5871	3.5605
	var	6.2670	6.5347	2.6140e+03	7.3790	7.0367	15.6203
SC1F	SR	0.0380	0.0902	−0.1253	−0.0510	−0.0034	0.6570
	TURN	0.3032	1.1648	0.8972	0.2964	0.2962	1.5858
	ASP	0.4690	1.3556	2.1771	0.5722	0.5863	2.8086

reports the OOS performance by using four quantities defined in Section 6.1. We set $\alpha =\beta =10$ and the sparsity level $s_1==\lceil 30\%n\rceil$ (on the DJIA, NASDAQ, SP500, and FF100 data sets) and $s_2==\lceil 30\%n\rceil$ (on the Russell2000 and Russell3000 data sets). We can observe that the RSA and RSQ portfolios achieve the smallest variances across all portfolio strategies, i.e., on average with $10.84{(\%)}^2$ and $11.09{(\%)}^2$, respectively. This means they are less volatile, i.e., less risky. SU, SC1F, and SCID have the highest variance on average, $995.73{(\%)}^2$, $442.81{(\%)}^2$ and $404.20{(\%)}^2$ in this setting. The variance of the remaining portfolio strategies is $11.94{(\%)}^2$ (L1), $11.98{(\%)}^2$ (EN), $16{(\%)}^2$ (L12), $27.78{(\%)}^2$ (SC), $28.71{(\%)}^2$ (EW), and $71.11{(\%)}^2$ (Box), respectively. In addition, we observe that the Sharpe ratios of the various portfolios on average are 12.34% (SC), 11.82% (EW), 11.77% (RSA), 11.70% (RSQ), 11.61% (L12), 11.21% (EN), 11.18% (L1), 10.58% (SCID), 10.09% (SC1F), 9.23% (Box), and 7.92% (SU). We see that the RSA and RSQ portfolios do not result in a significantly different OOS Sharpe ratio when compared with SC and EW; however, they are higher than the rest of the portfolio strategies.

As for the portfolio turnover, unsurprisingly, the EW portfolio strategy exhibit the lowest turnover of all portfolio strategies, amounting to 3.20%. The RSA and RSQ portfolio strategies have moderate levels of turnover on average, 13.49% and 13.20%. The highest average turnover is generated by the SU portfolio and, then, by the Box portfolio, amounting on average to 225.81% and 193.10%, meaning that they are very costly. The turnover of the remaining portfolio strategies range between 11.85% (L12), 25.76% (L2), 16.66% (L1), and 13.98% (EN), respectively. The high turnover of SU and Box was reflect in the enormous average short positions of over 283.52% and 333.52% on average across the six data sets. The second two highest average short positions are by SCID and SC1F, respectively, amounting to 164.29% and 132.81%. The average short positions of the SC and EW portfolios are on average approximately 0% across the six data sets. The average short positions of the RSQ and RSA portfolio strategies also tend to zero. Therefore, considering the moderate turnover and the average short positions, the proposed RSQ and RSA strategies represent a practically implementable method that outperform the portfolio strategies listed in Table 1.

6.4. Cumulative Return

In this subsection, we show the CR of several portfolio strategies. We use the FF100 data set. The sparsity level $s=\lceil 10\%n\rceil$. The parameter $\alpha =\beta =10$. According to the OOS performance, we choose RSQ, RSA, L12, L1, EN, EW, and SC to compare the CR.

Figure 4 shows the curves of the CR over the corresponding investment periods for the different portfolio strategies. Apparently, RSQ and RSA outperform the others with visible margins. However, RSA and RSQ do not produce significant differences. This result suggest that, compared with the other portfolios, the sparse portfolios RSA and RSQ grow more steadily together with a reduced volatility across most of the investment periods.

7. Conclusions

Portfolio selection has been a fertile area for robust optimization techniques. We proposed a robust and sparse portfolio selection optimization model by considering the perturbation in the asset return matrix and the parameter uncertainty in the expected asset return. We used the equivalence of robustness and regularization to deal with the perturbation in the asset return matrix. To deal with the uncertainty in the expected asset return, we considered two kinds of uncertainty sets and solved the worst-case scenario. We defined three types of stationary points of the penalty problem and then analyzed the relationship between these stationary points and local/global minimizers. Then, we designed the penalty alternating direction method to solve each problem. Although there is no theoretical guarantee for the equivalence between problems (5) and (6), as well as problems (8) and (18), we confirmed the effectiveness of our approach by comparing with nine widely studied portfolio models on seven real-world data sets. Extensive numerical experiments confirm that the portfolios we proposed have lower volatility, that is less risk. Moreover, our portfolio strategies can yield higher Sharpe ratios when the appropriate parameters are selected.

We note that the robust optimization (RO) mainly consider the uncertainty sets of parameters and thus it do not consider any distribution information of the data. This characteristic makes RO attractive, but at the same time, this method loses the comprehensive characterization of the data. Recently, distributed robust optimization (DRO) has attracted widespread attention and research. Although DRO takes into account the distribution information of the data, the cost paid is that it is difficult to solve. We will consider how to apply DRO to sparse portfolio problems, while considering the distribution information of financial data and improving the sparsity. The most direct extension is the distributed robust portfolio optimization with the ${\ell }_0$ norm constraint, which is a worthwhile and challenging issue.