ﬁfi ﬀff

Some connections between higher moments portfolio optimization methods
Farshad Noravesh111Email<EMAIL_ADDRESS>Kristiaan Kerstens
222Email<EMAIL_ADDRESS>
In this paper, different approaches to portfolio optimization having higher
moments such as skewness and kurtosis are classified so that the reader can
observe different paradigms and approaches in this field of research which is
essential for practitioners in Hedge Funds in particular. Several methods
based on different paradigms such as utility approach and multi-objective
optimization are reviewed and the advantage and disadvantageous of these ideas
are explained. Keywords: multi-objective optimization, portfolio optimization,
scalarization, utility

##  1 Introduction

One can see the huge articles on higher moments of portfolio optimization as a
spectrum of assumptions which varies from the most general and weakest
assumptions to the most restricted and strong assumptions. Classification of
these articles can be done in different contexts. An important context is
whether we use parametric modeling or nonparametric modeling. It is
interesting that these two paradigms have also been discussed in other fields
of research such as machine learning where kernel methods such as kernel
support vector machine(kernel-SVM) or nonlinear kernel dimensionality
reduction is a nonparametric approach while methods like variational auto-
encoder(VAE) is an example of parametric approaches. These two different
paradigms has its root in statistics where we have parametric and
nonparametric approach to solve statistical problems. An abstract way of
explaining these two paradigms is well understood of we consider parametric
approach as a way to make a finite dimensional approximation of probability
distributions of returns, while in the nonparametric approach, no explicit
form of distribution is assumed and our ignorance on reality avoids us to make
strong assumptions of reality. In reality, space of distributions on asset
returns are extremely complex. Thus, using parametric distributions builds a
finite dimensional structure on the space of return distributions and
therefore on the moments of the returns which are used in the objective
function of portfolio optimization. The importance of utility function is well
understood specially if a taylor series expansion of utility around the
expected return is written to realize the connection of utility with higher
moments such as skewness and kurtosis. Thus,

$u(R)=u(E(R))+u^{\prime}[E(R)][R-E(R)]+\frac{u"[E(R)]}{{2!}}\sigma^{2}_{R}+\frac{u^{\prime\prime\prime}[E[R]]}{3!}M^{3}_{R}$
(1.1)

where $\sigma^{2}_{R}$ is the variance of expected return, and is the third
moment about the mean of expected return which is skewness. The following
utility functions are used in the literature:

$\begin{split}u(R)&=ln(\lambda R)\\\ u(R)&=\sqrt{\lambda R}\\\
u(R)&=-e^{-\lambda R}\end{split}$ (1.2)

However, the third utility function is used more relatively for its nice
properties. Taking the expectation of utility after Taylor expansion would
produce:

$\begin{split}E(u(R))&=-e^{\lambda
E(R)}(1+\frac{\lambda}{1!}E(R-E(R))-\frac{\lambda^{2}}{2!}(R-E(R))^{2}+\frac{\lambda^{3}}{31!}E(R-E(R))^{3}\\\
&-\frac{\lambda^{4}}{4!}(R-E(R))^{4}+O(R^{5}))\end{split}$ (1.3)

Thus, the following optimization problem is produced:

$min_{\omega}-\omega^{\prime}\mu+\lambda_{1}\omega^{\prime}\Sigma\omega-\lambda_{2}\omega^{\prime}M_{3}(\omega\otimes\omega)+\lambda_{3}\omega^{\prime}M_{4}(\omega\otimes\omega\otimes\omega)$
(1.4)

where $\omega\otimes\omega$ stands for Kronecker product and,
$\lambda_{1}=\frac{\lambda^{2}}{2!}$ , $\lambda_{2}=\frac{\lambda^{3}}{3!}$ ,
$\lambda_{3}=\frac{\lambda^{4}}{4!}$ (Glawischnig & Seidl, 2013) has used an
iterative approach to solve this objective function without using polynomial
goal programming(PGP). They start with $\lambda$ equal to twenty and then
reduce it in each step until two to avoid short selling of more than one
hundred percent. The algorithm is multi-start to get robust results. The trick
that (Glawischnig & Seidl, 2013) has used is that in the first iteration the
skewness and kurtosis is calculated and then fixed to make the full
optimization problem a valid quadratic optimization to be solved by quadratic
solvers, but this iterative approach that changes the nature of the
optimization problem and neglects the third order polynomials and 4th order
polynomials could result a big difference. (Glawischnig & Seidl,
2013),(Jondeau & Rockinger, 2006) almost use similar ideas and their approach
has an advantage over PGP since the weights of the optimization are naturally
formed by the taylor approximation of utility rather than arbitrary selection
that is used in PGP.

##  2 Parametric and nonparametric modeling

The motivation comes from many different intuitions such as constraining the
structure of return distributions of financial securities and approximate the
complex nature of asset returns. (Adcock, 2014) uses coherent multivariate
probability distribution assumption for asset returns to be able to use
stein’s lemma to create the mean-variance-skewness efficient hyper-surface. So
the assumption is that the asset returns are defined by the convolution of a
multivariate elliptically symmetric distribution and a multivariate
distribution of non-negative random variables so that the efficient portfolios
could be computed using quadratic programming on the efficient surface. Many
researchers are amazed by the Multivariate skew normal (MSN) which was first
introduced in (Azzalini, 1996) and was applied in finance in (Adcock & Shutes,
1999) is the center of many portfolio selection methods. One of the
motivations of using MSN is the simplicity of the maximization of utility,
specially if utility is an exponential function of return as explained in
(Adcock, 2005),(Landsman et al., 2019), the drawback of such an approach is
explained in (Adcock, 2005) when the preference parameter is too high or very
small. But the bigger disadvantage of such an approach is simply shaping the
mean-variance-skewness efficient surface by choosing an exponential utility
function. Another drawback is that the generalization of this approach for
kurtosis and higher moments is not straightforward. A generalization of MSN
could be seen in (SAHU et al., 2003) where an analytic forms of densities are
obtained and there distributional properties are studied. These parametric
modeling approaches to model distributions for returns are very diverse and
many researchers have noticed that such as (Jondeaua & Rockinger, 2003) that
uses a generalized Student-t distribution or (Mencia & Sentan, 2009) which
uses location scale mixture of normals and using maximum likelihood to infer
the parameters. The log-normal distribution has also been used in the
literature to model asset returns, but its skewness is a function of the mean
and the variance, not a separate skewness parameter. An interesting example
which combines parametric modeling(log normal) of portfolio distribution with
goal programming is (Chang et al., 2008a). (Glawischnig & Seidl, 2013),Jondeau
& Rockinger (2006) are examples of nonparametric approach to portfolio
optimization in higher moments since no distribution is assumed. The problem
of skewness term in the optimization is that it makes the optimization a non-
convex problem and therefore not tractable. This is the motivation of (Konno
et al., 1993) that without any assumption on the form of probability
distribution of returns tries to approximate the third order term due to
skewness by a piecewise linear approximation, however this approximation
should be tested on more experimental data to see to what extent the
approximation is valid. (Konno et al., 1998) uses a similar approach and
resolves the third order nonlinearity of skewness by representing it with a
difference of two convex functions and then using branch and bound to solve
the mean-variance-skewness portfolio optimization problem. A more detailed and
visual derivation of (Konno et al., 1998) can be seen in (Konno & Suzuki,
1995),(Konno & Yamamoto, 2005). The problem with ideas in (Konno et al.,
1993),(Konno et al., 1998),(Konno & Suzuki, 1995),(Konno & Yamamoto, 2005) is
that it can not be generalized to higher moments easily and the order of
complexity will be high. Another example for nonparametric higher moment
portfolio optimization is based on the concept of shortage function and the
geometric representation of mean-variance-skewness portfolio is illustrated in
(Kerstens et al., 2011)

##  3 A priori, interactive and posteriori

The literature on higher moment portfolio optimization can be classified in a
different context. In this context, it is important how the preferences are
presented in the optimization process and therefore three different paradigms
are discussed in the literature namely a priori, interactive and posteriori. A
priori methods such as goal programming and utility function method are used
when the preference of investor are known beforehand. In goal programming,
first all objective optimization problems are solved regardless of other
objectives. Then a final scalar objective optimization is solved that has some
weights that need to be fixed. Even if all possible weights are checked, still
some Pareto efficient solutions may be missing. Examples of using goal
programming for higher moments portfolio optimization are (Aksarayli & Pala,
2018),(Bergh & van Rensburg, 2008). Although it is a simple algorithm, using
appropriate weights is a debate. It can also produce solutions that are not
Pareto efficient, and some algorithm is needed to project them back to the
Pareto efficient solutions. Goal programming has many variants as explained in
(Aouni et al., 2014),(Tamiz et al., 2013). 1- Lexicographic 2- weighted, see
for example (Chang, 2011) 3- polynomial, see for example (Chang et al.,
2008b),(Chunhachinda et al., 1997),(Davies et al., 2009),(Lai, 1991),(Mhiri &
Prigent, 2010),(Proelss & Schweizer, 2014) and also (Ghahtarani & Najafi,
2013) if robustness with respect to some coefficients of the optimization is a
concern as well. 4- stochastic 5- fuzzy Another paradigm is to use posteriori
methods, where the target is finding all Pareto efficient frontier. Methods
such as linear weighting method, weighted geometric mean and Normal Boundary
Intersection (NBI)(Audet et al., 2008), Modified Normal Boundary Intersection
(MNBI), Normal Constraint, Multiple-objective Branch-and-Bound,epsilon-
constraints method and Pascoletti Serafini scalarization. Finally, there is
interactive paradigm where the investor iteratively solves the optimization
and gets feedback from solutions to find the Pareto optimal solutions. An
example of using epsilon-constraints for portfolio optimization can be seen in
(Xidonas, Mavrotas & Psarras, 2010),(Xidonas, Mavrotasand & Psarras,
2010),(Xidonas et al., 2011) which also includes an interactive filtering
process to consider investor preferences. Methods like NBI are computationally
expensive since for each iteration, an optimization problem should be solved
but on the other hand has a geometrical intuition which can also be used for
other methods as is explained in (Ghane-Kanafi & Khorram, 2015)

##  4 Pascoletti Serafini Scalarization (SP)

SP has two parameters namely a and r, which are chosen from. The method is
originally from (Pascoletti & Serafini, 1984) but is generalized in
(Eichfelder, 2008). The novelty of their approach is on conditions to bound
the parameter a on a restricted hyperplane, but it can only be used in two
dimensions and therefore not applicable to our four dimensional objective
function that covers skewness and kurtosis as well. (Khorram et al., 2014)
attempts to find the restricted set for parameter a but only the trivial point
of zero for parameter a is considered and also the simple EP cone is
considered.

Figure 1: equivalence of multiobjective optimization algorithms

###  4.1 Shortage Function (SF)

(Briec et al., 2004) introduced the shortage function as a portfolio
performance measure in the traditional mean-variance portfolio framework. The
original shortage function is computed by solving the following problem:

$\begin{split}max\ \delta\\\ E[R(y^{k})]+\delta g_{E}&\leq E[R(x)],\\\
Var[R(y^{k})]-\delta g_{V}&\geq Var[R(x)],\\\ Sk[R(y^{k})]+\delta g_{S}&\leq
Sk[R(x)],\\\ \sum x_{i}&=1\ ,x_{i}\geq 0,i=1,...,n\end{split}$ (4.1)

where g is the direction vector and the effect of choosing g is well described
in (Kerstens et al., 2012). A good explanation of how shortage function can be
useful for mean-variance-skewness portfolio framework is explained in (Briec
et al., 2007). The connection between shortage function and NBI method is
explained in the present paper, however the connection between shortage
function and polynomial goal programming (PGP) is well described in (Briec et
al., 2013). Now, the following special case of shortage function is defined
which is called modified shortage function (MSF) and is useful to show its
connection with NBI method. subsection

######  Definition 1.

Modified shortage function(MSF) is defined as the solution of the following
optimization problem

$\begin{split}max\ \delta\\\ E[R(y^{k})]+\delta g_{E}&=E[R(x)],\\\
Var[R(y^{k})]-\delta g_{V}&=Var[R(x)],\\\ Sk[R(y^{k})]+\delta
g_{S}&=Sk[R(x)],\\\ \sum x_{i}&=1\ ,x_{i}\geq 0,i=1,...,n\end{split}$ (4.2)

###  2 NBI and SP

In this section, the equivalence of NBI and MSF are proved. One of the most
important scalarization technique is normal boundary intersection method (NBI)
which is as follows

$\begin{split}max\ s\\\ \Phi\beta+s\bar{n}&=f(x)-f^{*}\\\ s\in R,&\
x\in\Omega\end{split}$ (4.3)

Different optimization problems in (4.3) for different $\beta\in R^{m}_{+}$,
having $\sum_{i=1}^{m}\beta_{i}=1$ are solved. Here $f^{*}$ denotes the so-
called ideal point and the matrix $\Phi\in R^{m\times m}$ consists of the
columns $f(x^{i})-f^{*}\ ,i=1,\ldots,m$ and the vector $\bar{n}$ is defined as
normal unit vector to the hyperplane directing to the negative orthant. The
problem of this method is that not all minimal points can be found as a
solution of NBI. The following lemma shows the direct connection between NBI
and modified SP.

######  Lemma 2.1.

(Eichfelder, 2008) A point ($\bar{s},\bar{x}$) is a maximal solution of NBI
with $\beta\in R^{m}$ , $sum_{i=1}^{m}\beta_{i}=1$, if and only if
$(-\bar{s},\bar{x})$ is a minimal solution of $\overline{SP}(a,r)$ with
$a=f^{*}+\Phi\beta$ and $r=-\bar{n}$

###  3 NBI and MSF

Another connection is between NBI and modified shortage function as is proved
in the next proposition:

######  Proposition 3.1.

Modifed shortage function scalarization is equivalent to NBI by the following
substitutions:

$\begin{split}s&=\delta\\\ f_{1}(x)&=E[R(x)]\\\ f_{2}(x)&=V[R(x)]\\\
f_{3}(x)&=Sk[R(x)]\\\ \Phi\beta+f^{*}&=c\\\ \bar{n}&=g\end{split}$ (4.4)

where $c=\begin{bmatrix}E[R(y^{k})]\\\ Var[R(y^{k})]\\\
Sk[R(y^{k})]\end{bmatrix}$

###### Proof.

proof is straightforward by a simple substitution. ∎

###  4 SP and SF

######  Proposition 4.1.

Shortage function scalarization is equivalent to pascoletti serafini
scalarization by the following substitutions:

$\begin{split}\delta&=-t\\\
a_{1}&=E(R(y^{k^{\prime}}))=2E(R(x))-E(R(y^{k}))\\\
a_{2}&=V(R(y^{k^{\prime}}))=2V(R(x))-V(R(y^{k}))\\\
a_{3}&=E(R(y^{k^{\prime}}))=2S(R(x))-S(R(y^{k}))\\\ r_{1}&=g_{E}\\\
r_{2}&=g_{V}\\\ r_{3}&=g_{S}\end{split}$ (4.5)

###### Proof.

Since $max\ -t=min\ t$ and the reference point can be set as any point in
$R^{3}$, substituting the reference point(a) and the direction(r) in SP
constaint which is $a+tr-f(x)\in K$ and choosing the cone K to be trivial
$R^{3+}$ the proof is complete. ∎

###  5 NBI and goal programming

It is shown here that the popular goal programming(PGP) method which is widely
used in portfolio optimization literature has close connection with NBI and is
defined as follows: subsection

######  Definition 1.

Polynomial goal programming (PGP) is defined as:

$\begin{split}PGP(\alpha,\beta)&=min\\{d_{1}^{\alpha}+d_{3}^{\beta};d_{1}=z_{1}^{*}-z_{1},z_{2}=1,d_{3}=z_{3}^{*}-z_{3}\\}\\\
z_{1}^{*}&=max\\{z_{1};z_{2}=1\\}\\\
z_{3}^{*}&=max\\{z_{3};z_{2}=1\\}\end{split}$ (4.6)

######  Proposition 5.1.

A solution to NBI problem is also a solution to PGP portfolio optimization
problem.

###### Proof.

Since $(x^{*},s^{*},\lambda^{*})$ is the solution of NBI it satisfies the
first order KKT conditions:

$\begin{split}\nabla_{X}F(x^{*})\lambda^{*}&=0\\\
1+\hat{n}\lambda^{*}&=0\end{split}$ (4.7)

where $\lambda\in R^{3}$ represents the multipliers corresponding to the 3
equality constraints namely return,variance and skewness constraint. On the
other hand, the first oder KKT condition for PGP can be decomposed for two
different sets of coordinates. Stationarity equations for the first set of
coordinates $(\omega_{1},\omega_{2},\omega_{3})$ results:

$\nabla_{\omega}Z\ \mu^{\star}=0$ (4.8)

Now stationarity with respect to second set of coordinates $(d_{1},d_{3})$
yields:

$\begin{split}\alpha d_{1}^{\alpha-1}+\mu^{(1)}&=0\\\ \mu^{(2)}&=0\\\ \beta
d_{3}^{\beta-1}+\mu^{(3)}&=0\end{split}$ (4.9)

Now expanding $1+\hat{n}\lambda^{star}=0$ in 4.7 generates:

$1+\hat{n}_{1}\lambda_{1}^{*}+\hat{n}_{2}\lambda_{2}^{*}+\hat{n}_{3}\lambda_{3}^{*}=0$
(4.10)

Simplifying 4.8 would produce:

$1+\alpha\frac{d_{1}^{\alpha-1}}{\mu^{(1)}}=0$ (4.11)

combining (4.11) and (4.10) and solving for $\alpha$ , $\beta$ results:

$\begin{split}\alpha&=\frac{\hat{n}_{1}\lambda_{1}^{*}+\hat{n}_{2}\lambda_{2}^{*}+\hat{n}_{3}\lambda_{3}^{*}}{d_{1}^{\alpha-1}}\\\
\beta&=-\frac{\mu{(3)}}{d_{3}^{\beta-1}}\end{split}$ (4.12)

Equivalently, $(\omega^{*},d^{*},\mu^{*})$ is the solution of PGP problem and
since F and X and in NBI problem are the same as Z and $\omega$ in PGP problem
and for any solution and lagrange multipliers there exists equivalent ones by
suitable substitutions for $\alpha$ and $\beta$ like (4.12) the proof is
complete. ∎

##  5 Methodology

There are some quality measures that can be used to figure out which algorithm
is better than the other. These measures are well explained in (Eichfelder,
2008) namely coverage error($\epsilon$), uniformity level($\delta$) and
cardinality(number of points)
These three measures are conflicting in nature and another multiobjective
optimization on top of the main optimization problem may be formed for a
rigorous analysis. Satisfying these type of measures are what most researchers
refer to as adative considerations. In another context, there are two
important aspects in any multiobjective optimzation namely accuracy and
diversity. The former forces the solutions to converge to pareto frontier
while the latter makes the efficient set equidistanced as much as possible.
There are many other measures in the literature such as hypervolume indicator
but implementing some of them makes the algorithm very slow and even
convergence of them are not proved.

Figure 2: multiobjective optimization methods

So far, mostly Pascoletti Serafini methods are discussed in the present paper
but there are other ideas in the literature as depicted in Figure 2. Set
oriented methods steer a set of solutions at each iterations such as
(Hernandez et al., 2013) for cell mapping method or (Dellnitz et al., 2005)
for subdivision algorithms. In the present paper, two adaptive algorithm for
higher moment multiobjective portfolio optimization are given. The first one
is an adaptive epsilon constraint method while the second one is not based on
scalarization and has its root in (Hillermeier, 2001) but is recently
developed by many researchers as in (Martín & Schutze, 2018) , (Schutze et
al., 2020) and also generalized for problems having inequality constraints in
(Beltran et al., 2020). Both of the proposed algorithms in this section are
based on KKT conditions but the approaches are slightly different. Both
methods are designed to produce equidistance pareto frontier points. The
equidistance parameter in the first algorithm is $\alpha$ while in the second
algorithm it is called $\tau$ to mimick the variables in the related
historical articles. So both methods are adaptive in the sense of equal
distance points on efficient frontier but no other considerations are taken
for the delta and cardinality quality measures since they are expected to
produce good results, otherwise they add to the complexity of the
algorithms.The first algorithm is based on Epsilon constraint which is a
special case of Pascoletti Serafini while the second algorithm is based on
continuation methods and the connections are well shown in Figure 2.

###  1 Adaptive Epsilon Constraint

Adaptive ECS is described as

$\begin{split}min&\ f_{k}(x)\\\
f_{i}(x)&\leq\epsilon_{i},i\in\\{1,\ldots,m\\}\setminus\\{k\\}\\\
x\in\Omega\end{split}$ (5.1)

The scalar optimization problem (5.1) can be formulated as

$\begin{split}min\ t\\\ \epsilon_{i}-f_{i}(x)&\geq 0\
i\in\\{1,\ldots,m\\}\setminus\\{k\\}\\\ t-f_{k}(x)&\geq 0\\\ g_{j}(x)&\geq 0\
j\in\\{1,\ldots,p\\}\\\ h_{l}(x)&=0\ l\in\\{1,\ldots,q\\}\\\ t\in R&,\ x\in
R^{n}\\\ \end{split}$ (5.2)

It is proved by Theorem 2.27 in (Eichfelder, 2008) how it is possible to
relate SP to epsilon-constraint method via lagrange multipliers. So using the
following substitutions for $a_{i}$ and r, if $\bar{x}$ is a minimal solution
of (5.2) . Thus, (5.2) is equivalent to SP(a,r) with the following
substitutions for $a_{i}$ and r :

$\begin{split}a_{i}&=\epsilon_{i}\ \forall
i\in\\{1,\ldots,m\\}\setminus\\{k\\}\\\ a_{k}&=0\\\ r&=e_{k}\end{split}$ (5.3)

The full algorithm is shown in algorithm 1.

Algorithm 1 Adaptive $\epsilon$ constraint method for mean-variance-skewness

Input: Choose the desired number $N_{1}$ of discretization points for the
range of the functions $f_{1}$(i.e. in direction $v^{1}=(1,0,0)^{T}$ and
$N_{2}$ for the range of function $f_{2}$ (i.e. in direction
$v^{2}=(0,1,0)^{T}$ )
Step 1: solve the optimization problems $min_{x\in\Omega}f_{i}x$ with minimal
solution $x^{min,i}$ and minimal value $f_{i}(x^{min,i})$ for i=1,2 as well as
$max_{x\in\Omega}f_{i}x$ with maximal solution $x^{max,i}$ and maximal value
$f_{i}(x^{max,i})=\epsilon_{i}^{max}$ for i=1,2
Step 2: Set $L_{i}=\frac{\epsilon_{i}^{max}-\epsilon_{i}^{min}}{N_{i}}$ and
solve the problem $P_{3}(\epsilon)$
for all parameters $\epsilon\in E$ with
$E:=\\{\epsilon=(\epsilon_{1},\epsilon_{2})\in
R^{2}|\epsilon_{i}=\epsilon_{i}^{min}+\frac{L_{i}}{2}+l_{i}.L_{i}$
for $l_{i}=0,...,N_{i-1},i=1,2\\}$
Determine the set $A^{E}=\\{(\epsilon,\bar{x},\bar{\mu})|\bar{x}\\}$ where
$\bar{x}$ is a minimal solution of $(P_{3}(\epsilon))$ with parameter
$\epsilon$ and lagrange multiplier $\bar{\mu}$ to the constraints
$f_{i}(x)\leq\epsilon_{i}$
$,i=1,2$ for $\epsilon\in E$
Step 3: Determine the set $D^{H^{0},f}:=\\{f(x)|\exists\epsilon\in
R^{2},\mu\in R^{2}_{+}with(\epsilon,x,\mu)\in A^{E}\\}$
Input: Choose $y\in D^{H^{0},f}$ with $y=f(x^{\epsilon})$ and
$(\epsilon,x^{\epsilon},\mu^{\epsilon})\in A^{E}$ if y is a sufficient good
solution, then stop. Else, if additional points in the neighborhood of y are
desired, give a distance $\alpha\in R$, $\alpha>0$ in the image space and the
number of desired new points $\bar{n}=(2k+1)^{2}-1$ (for a $k\in N$) and go to
Step 4.
Step 4: Set
$\epsilon^{ij}:=\epsilon+i.\frac{\alpha}{1+(\mu_{1}^{\epsilon})^{2}}\begin{bmatrix}1\\\
0\\\
\end{bmatrix}+j.\frac{\alpha}{1+(\mu_{2}^{\epsilon})^{2}}\begin{bmatrix}0\\\
1\\\ \end{bmatrix}$
for all
$(i,j)\in\\{(i,j)\in Z^{2}|i,j\in\\{-k,\ldots,k\\},(i,j)\neq(0,0)\\}$

and solve problem $(P_{3}(\epsilon^{i,j}))$
if there exists a solution $x^{i,j}$ with Lagrange multiplier $\mu^{i,j}$,
then set $A^{E}:=A^{E}\cup\\{\epsilon^{i,j},x^{i,j},\mu^{i,j}\\}$.
Go to step 3.
Output: The set $D^{H_{0},f}$ is an approximation of the set of weakly
efficient points

Figure 3: Pareto Front using adaptive epsilon constraint method with 2500
points

The simulation results for adapative epsilon constraint method is illustrated
in Figure 3. As is explained in the previous section $\epsilon$ constraint
method is just a special case of SP and the following algorithm is used for
multiobjective portfolio optimization having three objectives namely return,
variance and skewness.

###  2 Adaptive Multi-start Pareto Tracer

There are two main ideas in this approach: 1-KKT condition in single objective
is generalized to multiobjective as is described in (Hillermeier, 2001). 2-A
predictor corrector idea which first predicts the next move in decision space
and then corrects it by a multiobjective gradient descent. The algorithm 2 is
a modification of (Martín & Schutze, 2018) by doing it in a multi start way
and shaping the objective space distribution and customizing it for portfolio
optimization of three objectives namely mean, variance and skewness. Consider
the multiobjective optimization problem defined bellow:

$\begin{split}min\ &F(x)\\\ s.t\ h(x)&=0\end{split}$ (5.4)

where F is a vector of objectives $F:R^{n}\rightarrow R^{k}$. F in 5.4 is
actually a three dimensional vector of mean,variance and skewness and decision
space has dimension n which refers to number of assets or factors in a
multifactor investment framework. h in 5.4 is the constraint that the sum of
allocations to different assets should be one. A predictor corrector method is
developed in (Hillermeier, 2001) by considering

$\tilde{F}(x,\alpha)=\begin{bmatrix}\sum_{i=1}^{k}\alpha_{i}\nabla f_{i}(x)\\\
\sum_{i=1}^{k}\alpha_{i}-1\end{bmatrix}=0$ (5.5)

The set of KKT points of 5.4 is contained in the null set of $\tilde{F}$ which
is the idea behind many continuation methods along $\tilde{F}^{-1}(0)$ as
written in 5.5. So a simple representation for the tangent vectors to Pareto
set can be written as

$\tilde{F}(x,\alpha)\begin{bmatrix}\nu\\\
\mu\end{bmatrix}=\begin{bmatrix}\sum_{i=1}^{k}\alpha_{i}\nabla^{2}f_{i}(x)&\nabla
f_{1}(x)&\ldots&\nabla f_{k}(x)\\\
\sum_{i=1}^{k}\alpha_{i}-1&1&\ldots&1\end{bmatrix}\begin{bmatrix}\nu\\\
\mu\end{bmatrix}=\begin{bmatrix}0\\\ 0\end{bmatrix}$ (5.6)

Thus the vector in 5.6 can be expressed as

$v_{\mu}=-W^{-1}_{\alpha}J^{T}\mu$ (5.7)

where $W_{\alpha}$ in (5.7) is defined as

$W_{\alpha}:=\sum_{i=1}^{k}\alpha_{i}\nabla^{2}f_{i}(x)\in R^{n\times m}$
(5.8)

and J is defined by

$J=J(x)=\begin{bmatrix}\nabla f_{1}(x)^{T}\\\ \vdots\\\ \nabla f_{k}(x)^{T}\in
R^{k\times n}\end{bmatrix}$ (5.9)

Now d has a special meaning that expresses the first order approximated
movement in objective space for infinitesimal step sizes and is defined below

$d:=J\nu$ (5.10)

Using (5.7) and the definition of d in (5.10) makes

$J\nu_{\mu_{d}}=-JW^{-1}_{\alpha}J^{T}\mu=d$ (5.11)

Now any possibility for selecting d will generate a different distribution on
pareto front. Since the resulting $\nu_{mu}$ are tangents to Pareto set and
the aim is making even spread of points on Pareto front, directions d should
be selected such that it makes orthonormal basis of tangent space to Pareto
front at F(x). One of the natural ways to do this is utilizing QR
factorization of $JW^{-1}J^{T}$ and by selecting

$d_{i}:=q_{i+1}\ i=1,\ldots,k-1$ (5.12)

$\mu$ can be solved. Now it is possible to obtain the predictor $p:=x+t\nu$ to
make an evenly distributed set of solutions along the Pareto front the
following approximation for t can be chosen

$\begin{split}\lVert F(x_{i})-F(x_{i+1})&\approx\tau\rVert\\\
t=\frac{\tau}{\lVert J\nu_{\mu}\rVert}\end{split}$ (5.13)

The predictor part explained so far cares about going through the tangent
direction in Pareto set to create an equidistanced set of Pareto front points.
The next part is the corrector which cares about convergence issues and is
explained in (Fliege et al., 2009)and the present paper implements it for the
corrector part of the algorithm, but there are other methods that can be used
for this step such as (Povalej, 2014) which approximates the second derivative
matrices instead of evaluating them although both methods have superlinear
rate of convergence. The corrector part implemented in the present paper is
based in minimization of the following optimization problem.

$\begin{split}min\ g(t,s)&=t\\\ s.t\ \nabla
F_{j}(x)^{T}s+\frac{1}{2}s^{T}\nabla^{2}F_{j}(x)s-t&\leq 0\ (t,s)\in R\times
R^{n}\end{split}$ (5.14)

Thus, the full predictor-corrector algorithm is shown in Algorithm 2.

Algorithm 2 Adaptive Multi-start Pareto Tracer for mean-variance-skewness

Repeat Predictor-corrector loop:
Input : create some bundles to prepare for the multistart algorithm Predictor
part:
Step 1: calculate $\mu$ by solving (5.11)
Step 2: calculate direction $v_{mu}$ from (5.7)
Step 3: update the predictor position in decision space by $p:=x+t\nu_{\mu}$
Corrector part:
Step 4: solving (5.14) and update x

Figure 4: adaptive multistart Pareto Tracer

Simulation results are shown in Figure 4. Since the current algorithm has many
parameters to tune and to get all parts of pareto front in a faster way, an
evolutionary algorithm could be combined with the current algorithm to make a
hybrid algorithm.

##  6 Conclusions

The paradigms in higher moment multiobjective optimization are critically
reviewed and the connection between some of them are explained in the present
paper. It has been proved that shortage function method can be seen as a
Pascoletti Serafini scalarization. Finally, two algorithms for portfolio
optimization are suggested. The first one is based on scalarization paradigm
and is called adaptive epsilon constraint method while the second one is a
type of continuation method and is called adaptive multistart Pareto Tracer
which bundles different local solutions to provide a global Pareto Front by
both exploration and exploitation.

##  7 Future Works

The first suggested algorithm can be modified by handling variable structure
ordering as is explained in (Eichfelder, 2014) and (Eichfelder, 2012). So
instead of a fixed cone K for ordering, a variable cone is considered and each
point could have a different ordering corresponding to a different cone. The
second suggested algorithm can be further developed by hybridizing it with
evolutionary algorithms such as genetic algorithm to make it faster. Another
line of research which could be theoretically interesting is to see if the
continuation method suggested in the second algorithm could be seen as
generalization of Pascoletti Serafini ,since the continuation framework is
directly working on KKT conditions for a multiobjective optimization problem
while scalarization takes advantage of KKT condition of a mono objective
problem.

## References

  * (1)
  * Adcock (2005) Adcock, C. (2005), ‘Exploiting skewness to build an optimal hedge fund with a currency overlay’, The European Journal of Finance 11, 445–462.
  * Adcock (2014) Adcock, C. (2014), ‘Mean-variance-skewness efficient surfaces’, European Journal of operational research 234, 392–401.
  * Adcock & Shutes (1999) Adcock, C. & Shutes, K. (1999), ‘Portfolio selection based on the multivariate skew normal distribution’, Financial Modeling 11, 167–177.
  * Aksarayli & Pala (2018) Aksarayli, M. & Pala, O. (2018), ‘A polynomial goal programming model for portfolio optimization based on entropy and higher moments’, Expert systems with applications 94, 185–192.
  * Aouni et al. (2014) Aouni, B., Colapinto, C. & Torre, D. L. (2014), ‘Financial portfolio management through the goal programming model: Current state-of-the-art’, European Journal of Operational Research 234, 536–545.
  * Audet et al. (2008) Audet, C., Savard, G. & Zghal, W. (2008), ‘Multiobjective optimization through a series of single objective formulations’, Journal of optimization 19, 188–210.
  * Azzalini (1996) Azzalini, D. (1996), ‘The multivariate skew normal distribution’, Biometrika 83, 715–726.
  * Beltran et al. (2020) Beltran, F., Cuate, O. & Schutze, O. (2020), ‘The pareto tracer for general inequality constrained multi-objective optimization problems’, Mathematical and Computational Applications 25(4).
  * Bergh & van Rensburg (2008) Bergh, G. & van Rensburg, P. (2008), ‘Hedge funds and higher moment portfolio selection’, Journal of Derivatives and Hedge Funds 14, 269–297.
  * Briec et al. (2007) Briec, W., Kerstens, K. & Jokung, O. (2007), ‘Mean variance skewness portfolio performance gauging: A general shortage function and its dual approach’, Management science 53, 135–149.
  * Briec et al. (2004) Briec, W., Kerstens, K. & Lesourd, J. (2004), ‘Single period markowitz portfolio selection, performace gauging and duality, a variation on luenberger shortage function’, Journal of optimization theory and applications 120, 1–27.
  * Briec et al. (2013) Briec, W., Kerstens, K. & Woestyne, I. (2013), ‘Portfolio selection with skewness: A comparison of methods and a generalized one fund result’, European journal of operational research 230, 412–421.
  * Chang et al. (2008a) Chang, C., DuPoyet, B. & Prakash, A. (2008a), ‘Optimum allocation of weights to assets in a portfolio: the case of nominal annualization versus effective annualization of returns’, Applied Financial Economics 18, 1635–1646.
  * Chang et al. (2008b) Chang, C.-H., Dupoyet, B. & Prakash, A. (2008b), ‘Effect of intervalling and skewness on portfolio selection in developed and developing markets’, Applied Financial Economics 18, 1697–1707.
  * Chang (2011) Chang, C.-T. (2011), ‘Multi-choice goal programming with utility functions’, European Journal of Operational Research 215, 439–445.
  * Chunhachinda et al. (1997) Chunhachinda, P., Dandapani, K. & Shahid Hamid, A. P. (1997), ‘Portfolio selection and skewness: Evidence from international stock markets’, Journal of banking and finance 21, 143–167.
  * Davies et al. (2009) Davies, R., Kat, H. & Lu, S. (2009), ‘Fund of hedge funds portfolio selection a multiple-objective approach’, Journal of Derivatives and Hedge Funds 15, 91–115.
  * Dellnitz et al. (2005) Dellnitz, M., Schutze, O. & Hestermeyer, T. (2005), ‘Covering pareto sets by multilevel evolutionary subdivision techniques’, Journal of Optimization Theory and Applications 124, 113–136.
  * Eichfelder (2008) Eichfelder, G. (2008), Adaptive Scalarization Methods in Multiobjective Optimization, springer.
  * Eichfelder (2012) Eichfelder, G. (2012), ‘Numerical procedures in multiobjective optimization with variable ordering structures’, Journal of Optimization Theory and Applications 42, 1–26.
  * Eichfelder (2014) Eichfelder, G. (2014), Variable Ordering Structures in Vector Optimization, springer.
  * Fliege et al. (2009) Fliege, J., Drummond, L. M. G. & Svaiter, B. F. (2009), ‘Newton’s method for multiobjective optimization’, SIAM Journal on Optimization 20, 602–626.
  * Ghahtarani & Najafi (2013) Ghahtarani, A. & Najafi, A. (2013), ‘Robust goal programming for multi objective portfolio selection problem’, Economic Modelling 33, 588–592.
  * Ghane-Kanafi & Khorram (2015) Ghane-Kanafi, A. & Khorram, E. (2015), ‘A new scalarization method for finding the efficient frontier in non-convex multi-objective problems’, applied mathematical modelling 39, 7483–7498.
  * Glawischnig & Seidl (2013) Glawischnig, M. & Seidl, I. (2013), ‘Portfolio optimization with serially correlated, skewed and fat tailed index returns’, CEJOR 21, 153–176.
  * Hernandez et al. (2013) Hernandez, C., Naranjani, Y., Sardahi, Y., Liang, W., Schutze, O. & Sun, J.-Q. (2013), ‘Simple cell mapping method for multi-objective optimal feedback control design’, International Journal of Dynamics and Control 1, 231–238.
  * Hillermeier (2001) Hillermeier, C. (2001), Nonlinear multiobjective optimization: A generalized homotopy approach, Vol. 25, International Series on Numerical Mathematics. Basel: Birkhauser.
  * Jondeau & Rockinger (2006) Jondeau, E. & Rockinger, M. (2006), ‘Optimal portfolio allocation under higher moments’, European Financial Management 12, 29–55.
  * Jondeaua & Rockinger (2003) Jondeaua, E. & Rockinger, M. (2003), ‘Conditional volatility, skewness, and kurtosis: existence, persistence, and comovements’, Journal of Economic Dynamics and Control 27, 1699–1737.
  * Kerstens et al. (2011) Kerstens, K., Mounir, A. & Woestyne, I. (2011), ‘Geometric representation of the mean–variance–skewness portfolio frontier based upon the shortage function’, European journal of the operational research 210, 81–94.
  * Kerstens et al. (2012) Kerstens, K., Mounir, A. & Woestyne, I. (2012), ‘Benchmarking mean-variance portfolios using a shortage function: The choice of direction vector affects rankings’, Journal of the operational research society 63, 1199–1212.
  * Khorram et al. (2014) Khorram, E., Khaledian, K. & Khaledyan, M. (2014), ‘A numerical method for constructing the pareto front of multi-objective optimization problems’, Journal of Computational and Applied Mathematics 261, 158–171.
  * Konno et al. (1993) Konno, H., Shirakawa, H. & Yamazaki, H. (1993), ‘A mean-absolute deviation-skewness portfolio optimization model’, Annals of operational research 45, 205–220.
  * Konno & Suzuki (1995) Konno, H. & Suzuki, K. (1995), ‘A mean-variance-skewness portfolio optimization model’, Journal of operations research 38, 173–187.
  * Konno et al. (1998) Konno, H., Suzukia, T. & Kobayashi, D. (1998), ‘A branch and bound algorithm for solving mean-risk-skewness portfolio models’, Optimization methods and software 10, 297–317.
  * Konno & Yamamoto (2005) Konno, H. & Yamamoto, R. (2005), ‘A mean-variance-skewness model: Algorithm and applications’, International Journal of Theoretical and Applied Finance 8, 409–423.
  * Lai (1991) Lai, T.-Y. (1991), ‘Portfolio selection with skewness’, Review of quantitative finance and accounting 1, 293–305.
  * Landsman et al. (2019) Landsman, Z., Makov, U. & Shushi, T. (2019), ‘Analytic solution to the portfolio optimization problem in a mean-variance-skewness model’, The European Journal of Finance 1, 165–178.
  * Martín & Schutze (2018) Martín, A. & Schutze, O. (2018), ‘Pareto tracer: a predictor–corrector method for multi-objective optimization problems’, Engineering Optimization 50, 516–536.
  * Mencia & Sentan (2009) Mencia, J. & Sentan, E. (2009), ‘Multivariate location–scale mixtures of normals and mean–variance–skewness portfolio allocation’, Journal of Econometrics 153, 105–121.
  * Mhiri & Prigent (2010) Mhiri, M. & Prigent, J. (2010), ‘International portfolio optimization with higher moments’, International journal of economics and finance 2, 157–169.
  * Pascoletti & Serafini (1984) Pascoletti, A. & Serafini, P. (1984), ‘Scalarizing vector optimization problems’, Optimization Theory and Applications 42, 499–524.
  * Povalej (2014) Povalej, Z. (2014), ‘Quasi-newton’s method for multiobjective optimization’, Journal of Computational and Applied Mathematics 255, 765–777.
  * Proelss & Schweizer (2014) Proelss, J. & Schweizer, D. (2014), ‘Polynomial goal programming and the implicit higher moment preferences of us institutional investors in hedge funds’, Financial Markets and Portfolio Management 28, 1–28.
  * SAHU et al. (2003) SAHU, S., DEY, D. & BRANCO, M. (2003), ‘A new class of multivariate skew distributions with applications to bayesian regression models’, The Canadian Journal of Statistics 31, 129–150.
  * Schutze et al. (2020) Schutze, O., Cuate, O., Martín, A., Peitz, S. & Dellnitz, M. (2020), ‘Pareto explorer: a global/local exploration tool for many-objective optimization problems’, Engineering Optimization 52, 832–855.
  * Tamiz et al. (2013) Tamiz, M., Azmi, R. & Jones, D. (2013), ‘On selecting portfolio of international mutual funds using goal programming with extended factors’, European Journal of Operational Research 226, 560–576.
  * Xidonas, Mavrotas & Psarras (2010) Xidonas, P., Mavrotas, G. & Psarras, J. (2010), ‘Equity portfolio construction and selection using multiobjective mathematical programming’, J Glob Optim 47, 185–209.
  * Xidonas et al. (2011) Xidonas, P., Mavrotas, G. & Zopounidis, C. (2011), ‘Ipssis: An integrated multicriteria decision support system for equity portfolio construction and selection’, European Journal of Operational Research 210, 398–409.
  * Xidonas, Mavrotasand & Psarras (2010) Xidonas, P., Mavrotasand, G. & Psarras, J. (2010), ‘Portfolio construction on the athens stock exchange: a multiobjective optimization approach’, Optimization 59, 1211–1229.