# Portfolio Optimization Rules beyond the Mean-Variance Approach

Maxime Markov and Vladimir Markov Corresponding Author:
<EMAIL_ADDRESS>
###### Abstract

In this paper, we revisit the relationship between investors’ utility
functions and portfolio allocation rules. We derive portfolio allocation rules
for asymmetric Laplace distributed $ALD(\mu,\sigma,\kappa)$ returns and
compare them with the mean-variance approach, which is based on Gaussian
returns. We reveal that in the limit of small $\frac{\mu}{\sigma}$, the
Markowitz contribution is accompanied by a skewness term. We also obtain the
allocation rules when the expected return is a random normal variable in an
average and worst-case scenarios, which allows us to take into account
uncertainty of the predicted returns. An optimal worst-case scenario solution
smoothly approximates between equal weights and minimum variance portfolio,
presenting an attractive convex alternative to the risk parity portfolio.
Utilizing a microscopic portfolio model with random drift and analytical
expression for the expected utility function with log-normal distributed
cross-sectional returns, we demonstrate the influence of model parameters on
portfolio construction. Finally, we address the issue of handling singular
covariance matrices by imposing block structure constraints on the precision
matrix directly. This comprehensive approach enhances allocation weight
stability, mitigates instabilities associated with the mean-variance approach,
and can prove valuable for both short-term traders and long-term investors.

##  1 Introduction

Investors and portfolio managers constantly strive to optimize their
investment strategies in order to maximize returns while minimizing risk. This
objective is formalized in the portfolio optimization rules, which facilitate
the selection of a set of stocks that yield the highest returns possible
considering the risk involved.

Since the introduction of the mean-variance approach in 1952 [1], the field
has rapidly evolved to include different risk and expected return models into
the framework. Despite this, simplistic rules, such as equal weight or equal
volatility, have demonstrated resilience in practice [2], while many advanced
formulas remain underutilized. The framework’s idealistic presumptions are the
issue. The mean-variance approach is based on the assumption that stock
returns are Gaussian, whereas real-world stock returns often exhibit fat tails
and skewness toward negative returns. Although expected returns display some
predictability over the long term because of their correlation with
fundamental data, mid- and short-term returns are inherently unpredictable.
Practically, the framework’s reliance on inverting the low-rank covariance
matrix introduces numerical instabilities to an already extensive list of
assumptions that are often invalid in a noisy and non-stationary world.
Critiques of modern portfolio theory span from popular authors [3] to
academics [4].

In this paper, we propose a set of models to enhance the quality of asset
allocation rules by incorporating uncertainty of the expected returns in
average and worst-case scenario and improving stability of covariance matrix
estimation by imposing block constraints on a precision matrix. First, we
examine historical returns of major stock indexes on daily, weekly, and
monthly time scales, demonstrating that the asymmetric Laplace distribution
(ALD) provides a reasonable fit. We derive empirical location ($\mu$), scale
($\sigma$), and skewness ($\kappa$) parameters as well as their scaling laws.
Second, we derive allocation rules for ALD returns and systematically compare
them to the mean-variance approach. We demonstrate that in the limit of small
$\frac{\mu}{\sigma}$, the Markowitz contribution is accompanied by a non-
negligible skewness term. Furthermore, the fat tails impose a cut-off on risk
by saturating allocation weights for stocks with high expected returns. Third,
we study the proposed models when the expected return $\mu$ is a Gaussian
random variable $\mu\sim N(\mu_{0},\sigma_{0})$, with the main objective to
find allocation rules when uncertainty $\sigma_{0}$ is much larger than
$\mu_{0}$. We derive allocation rules in an average scenario by
marginalization over $\mu$, and in the worst-case scenario. We demonstrate
that in the limit of noisy prediction of the expected return
$\frac{\mu_{0}}{\sigma_{0}}\to 0$, the allocation rule of the average scenario
resembles the mean-variance solution with regularization, where $\sigma_{0}$
acts as a regularizer (shrinkage). In the worst-case scenario, the optimal
solution smoothly approximates between equal weights and a minimum variance
portfolio, thus providing an attractive convex alternative to a risk parity
(RP) portfolio. Fourth, we derive allocation rules for log-normal distributed
returns, which serve as a base model for long-term cross-sectional returns
[5]. Finally, we investigate the problem of integrating the block structure of
investments into portfolio allocation rules. We show that block constraints
can be naturally incorporated into precision matrix estimation, resulting in
conditional independence between inter-block constituents. Block structure
mitigates the sampling error of large covariance matrices. It enables the
inversion of singular covariance or precision matrices if the maximum size of
a block is smaller than the number of observations in the sample.

We anticipate that the proposed framework and methods will prove beneficial
for practitioners in the realm of portfolio construction, potentially
improving or superseding the mean-variance and risk parity approaches.

##  2 Risk perception, utility theory, and portfolio allocation rules

Portfolio optimization task can be defined in many ways, as it is contingent
upon the diversity of risk perception and the distribution of investment
outcomes. In the standard economic approach, asset allocations are derived
from maximizing the expected value of the investor’s utility function. The
most commonly used utility models in finance are the Constant Absolute Risk-
Aversion (CARA) and Constant Relative Risk-Aversion (CRRA) models [6]. The
utility functions $U_{a}$ for CARA and $U_{r}$ for CRRA with respect to the
investment outcome (return) $x$ and with $a$ and $\gamma$ being the risk
aversion parameters can be expressed as follows:

$U_{a}(x,a)=\begin{cases}\frac{1-\exp(-ax)}{a},&\mbox{if }a\neq 0\\\
x,&\mbox{if
}a=0\end{cases},\,\,\,U_{r}(x,\gamma)=\begin{cases}\frac{x^{1-\gamma}-1}{1-\gamma},&\mbox{if
}\gamma>0,\gamma\neq 1\\\ \ln{x}&\mbox{if }\gamma=1\end{cases}$ (1)

To determine the optimal portfolio weights, we maximize the expected value of
the utility function $U(x,a)$ with respect to $w$:

$w=arg\max_{w}E_{P}[U(x,a)]=arg\max_{w}\int dx\,\,U(x,a)P(x)$ (2)

where $P(x)$ is a distribution of outcome $x$. Mean variance approach assumes
that returns $P(x)$ are Gaussian. In Section 3, we show that ALD provides a
better fit while capturing fatter tails and skewness of empirical returns.

##  3 Empirical distribution of stock returns on daily, weekly, and monthly
time scales

In this section, we examine the empirical return distribution $P(x)$ of 38
stock index time series over a period from January 1, 2002, to December 31,
2022, on daily ($D$), weekly ($W$), and monthly ($M$) time scales. We group
indexes according to their geographical location into several groups,
including the United States, Europe, Asia-Pacific (APAC), Japan, and BRIC
(Brazil, Russia, India, and China) countries. We also study 10 sectors of S&P
500 GICS Level 1 indexes to better understand the role of heterogeneity within
the S&P 500 index. We assume that index returns serve as a proxy for
individual stock returns over the chosen time frame and market.

All time series exhibit non-normal characteristics, specifically fat-tails and
skewness in return distribution. To account for these properties, we use the
ALD with the probability density function (PDF) given by [7]:

$p(x,\mu,\sigma,\kappa)=\frac{\sqrt{2}}{\sigma}\frac{\kappa}{1+\kappa^{2}}\exp{\left(-\sqrt{2}\frac{(x-\mu)}{\sigma}s\kappa^{s}\right)}$
(3)

where $s=sgn(x-\theta)$, $\mu$ is the location parameter, $\sigma$ is the
scale parameter and $\kappa$ is the skewness parameter. The case of zero
skewness corresponds to $\kappa=1$. The skewness toward negative returns
corresponds to $\kappa>1$. Note that there is a widely used classical ALD,
which differs from the ALD in Eq. 3 by $\sqrt{2}$ in the definition of scale
$\sigma$.

Figure 1: Probability density functions (PDF) of empirical daily (blue),
weekly (orange), and monthly (green) returns of S&P 500 and NASDAQ indexes.
Dashed lines show the mean posterior of the ALD fitted with the MCMC method.
Higher period distributions show greater asymmetry biased toward negative
returns.

The empirical PDF of the daily (blue), weekly (orange), and monthly (green)
returns of S&P500 (SPX) and NASDAQ (CCMP) indexes are shown in Figure 1. It
can be noticed that the daily returns have a sharp peak at zero. As the time-
frame increases, the distribution becomes broader and more skewed toward
negative returns. In order to quantify this effect, we fit the historical
return data with the ALD. The distribution parameters are estimated using
Markov Chain Monte Carlo (MCMC) method as implemented in PyMC3 [8]. The dashed
lines represent the mean posterior of the fits.

We repeat the same procedure for all indexes and their distribution parameters
$\mu$, $\sigma$, and $\kappa$ for the daily, weekly, and monthly returns are
summarized in Table 1. To highlight the effect of increasing skewness, we have
colored the columns with $\kappa$ values for each return.

Given the dependence of the distribution parameters for the indexes presented
in Table 1, we seek to extract the scaling laws for the averaged parameters.
To achieve this, we utilize a linear regression method on a log-log scale. By
analyzing the slope of a linear curve, we obtain the following scaling laws:

$\mu_{T}\sim T^{0.78},\,\,\sigma_{T}\sim T^{0.46}$ (4)

where $T$ represents the measured period in trading days. Both the location
and scale parameters exhibit a positive slope. In contrast, the behavior of
the skewness parameter $\kappa$ is idiosyncratic, and it may exhibit non-
monotonic behavior. An increase in $\kappa$ signifies the persistence of
negative returns across the investigated time scales. Therefore, it would be
prudent to select $\kappa$ based on the time scale that corresponds to the
anticipated portfolio holding period.

We observe a linear relation $\frac{\mu_{T}}{\sigma_{T}}=a*(\kappa_{T}-1)+b$,
with parameters $a_{D}=2.8$, $b_{D}=0.1$, $a_{W}=2.9$, $b_{W}=0.1$, and
$a_{M}=2.8$, $b_{M}=0.3$ for daily, weekly, and monthly data as summarized in
Table 1. The regression chart is shown in Appendix D.

Daily Weekly Monthly $\mu$ $\sigma$ $\kappa$ $\mu$ $\sigma$ $\kappa$ $\mu$
$\sigma$ $\kappa$ US SPX 0.10 0.57 high!17!low!90 1.04 0.41 1.19
high!35!low!90 1.08 1.48 2.17 high!57!low!90 1.13 CCMP 0.15 0.69
high!26!low!90 1.06 0.52 1.46 high!35!low!90 1.08 1.92 2.76 high!52!low!90
1.12 RIY 0.10 0.57 high!17!low!90 1.04 0.43 1.20 high!39!low!90 1.09 1.56 2.20
high!57!low!90 1.13 RTY 0.16 0.76 high!26!low!90 1.06 0.46 1.59 high!30!low!90
1.07 1.73 2.81 high!48!low!90 1.11 RAY 0.10 0.58 high!17!low!90 1.04 0.44 1.22
high!39!low!90 1.09 1.61 2.22 high!61!low!90 1.14 RLV 0.10 0.58 high!17!low!90
1.04 0.37 1.22 high!30!low!90 1.07 1.36 2.21 high!52!low!90 1.12 RLG 0.12 0.60
high!22!low!90 1.05 0.37 1.27 high!26!low!90 1.06 1.53 2.36 high!48!low!90
1.11 NBI 0.14 0.84 high!17!low!90 1.04 0.46 1.73 high!22!low!90 1.05 0.61 3.14
0.98 S&P 500 S5COND 0.14 0.67 high!26!low!90 1.06 0.29 1.46 high!13!low!90
1.03 0.82 2.68 high!4!low!90 1.01 S5CONS 0.07 0.44 high!17!low!90 1.04 0.41
0.95 high!43!low!90 1.10 1.09 1.79 high!39!low!90 1.09 S5ENRS 0.09 0.88
high!9!low!90 1.02 0.51 1.83 high!26!low!90 1.06 0.98 3.36 high!9!low!90 1.02
S5FINL 0.07 0.81 high!9!low!90 1.02 0.43 1.70 high!30!low!90 1.07 1.55 2.93
high!52!low!90 1.12 S5HLTH 0.10 0.54 high!17!low!90 1.04 0.44 1.17
high!39!low!90 1.09 1.07 2.08 high!26!low!90 1.06 S5INFT 0.15 0.74
high!22!low!90 1.05 0.57 1.58 high!35!low!90 1.08 1.97 2.99 high!48!low!90
1.11 S5MATR 0.14 0.75 high!22!low!90 1.05 0.58 1.56 high!43!low!90 1.10 0.94
2.93 high!13!low!90 1.03 S5RLST 0.14 0.80 high!22!low!90 1.05 0.45 1.64
high!30!low!90 1.07 1.21 2.87 high!30!low!90 1.07 S5TELS 0.09 0.67
high!17!low!90 1.04 -0.11 1.44 0.97 0.85 2.79 high!35!low!90 1.08 S5UTIL 0.20
0.58 high!48!low!90 1.11 0.30 1.26 high!22!low!90 1.05 1.56 2.17
high!65!low!90 1.15 S5INDU 0.13 0.65 high!22!low!90 1.05 0.34 1.42
high!22!low!90 1.05 1.09 2.59 high!26!low!90 1.06 Europe DAX 0.15 0.70
high!26!low!90 1.06 0.63 1.53 high!48!low!90 1.11 1.60 2.85 high!48!low!90
1.11 CAC 0.10 0.68 high!17!low!90 1.04 0.57 1.46 high!52!low!90 1.12 1.28 2.61
high!52!low!90 1.12 BEL20 0.08 0.61 high!17!low!90 1.04 0.68 1.33
high!74!low!90 1.17 1.69 2.35 high!83!low!90 1.19 UKX 0.09 0.56 high!22!low!90
1.05 0.42 1.15 high!48!low!90 1.11 1.14 1.99 high!61!low!90 1.14 IBEX 0.12
0.70 high!26!low!90 1.06 0.70 1.55 high!65!low!90 1.15 1.16 2.74
high!48!low!90 1.11 KFX 0.15 0.64 high!26!low!90 1.06 0.93 1.33 high!87!low!90
1.20 1.28 2.58 high!22!low!90 1.05 OMX 0.11 0.68 high!17!low!90 1.04 0.56 1.41
high!48!low!90 1.11 0.80 2.55 high!17!low!90 1.04 SMI 0.10 0.55 high!26!low!90
1.06 0.62 1.14 high!74!low!90 1.17 1.39 1.94 high!78!low!90 1.18 APAC AS51
0.123 0.50 high!35!low!90 1.08 0.49 1.04 high!61!low!90 1.14 1.54 1.96
high!78!low!90 1.18 Japan NKY 0.09 0.72 high!13!low!90 1.03 0.61 1.53
high!48!low!90 1.11 1.07 2.73 high!26!low!90 1.06 TPX 0.09 0.66 high!17!low!90
1.04 0.61 1.41 high!57!low!90 1.13 1.00 2.52 high!30!low!90 1.07 BRIC IBOV
0.13 0.88 high!13!low!90 1.03 0.64 1.89 high!30!low!90 1.07 0.53 3.46 0.96
RTSI 0.24 1.00 high!30!low!90 1.07 0.75 2.32 high!35!low!90 1.08 1.10 4.36
high!4!low!90 1.01 NIFTY 0.15 0.67 high!22!low!90 1.05 0.80 1.47
high!52!low!90 1.12 0.93 3.14 0.97 MXIN 0.18 0.66 high!26!low!90 1.06 0.87
1.47 high!61!low!90 1.14 0.91 3.09 0.97 SHCOMP 0.08 0.73 high!13!low!90 1.03
0.15 1.67 high!4!low!90 1.01 0.31 3.48 0.98 SHSZ300 0.11 0.81 high!13!low!90
1.03 0.40 1.85 high!17!low!90 1.04 0.57 3.77 0.98

Table 1: Location $\mu$, scale $\sigma$, and asymmetry $\kappa$ parameters of
ALD for daily, weekly and monthly returns (in percent) obtained from the MCMC
fits. Asymmetry values are highlighted with gradient colors varying from min
(yellow) to max (red) values.

##  4 Univariate portfolio allocation rules

In Section 3, we showed that the ALD, which has both fat tails and skewness,
provides a good fit for the empirical data. Given the distribution $P(x)$ and
the utility function $U(x,a)$, we can now determine the optimal portfolio
weights.

###  4.1 Univariate Gaussian returns

To develop intuition, we first discuss the mean-variance framework, which
assumes Gaussian returns. Given that the portfolio return $x$ follows a normal
distribution $x\sim N(\mu,\sigma)$, the expected utility function
$E_{N}[U_{a}(x,a)]$ can be expressed as:

$E_{N}[U_{a}(x,a)]=\int
dx\,\,U_{a}(x,a)N(\mu,\sigma)=\frac{1-e^{\frac{1}{2}a^{2}\sigma^{2}-a\mu}}{a}$
(5)

where we examine the CARA utility function $U_{a}(x,a)$ from Eq. 1. According
to Eq. 2, the optimal portfolio maximizes the expected utility function in Eq.
5. This is equivalent to maximizing the expression
$\mu-\frac{1}{2}a\sigma^{2}$.

Suppose there is one dollar to invest for a period of $t$ days. There are two
possible investment choices: a risky asset with a random return $r_{t}$ that
follows a probability distribution $N(\mu,\sigma)$, and a riskless asset with
a fixed return $r_{0}$. Let the portfolio allocation weight assigned to the
risky asset be $w$ and the weight assigned to the riskless asset be $(1-w)$.
Then, the investment result $x$ can be expressed as:

$x=w(1+r_{t})+(1-w)(1+r_{0})=1+r_{0}+w(r_{t}-r_{0})$ (6)

This expression follows a normal distribution with mean
$\mu_{x}=1+r_{0}+w(\mu-r_{0})$ and standard deviation $\sigma_{x}=w\sigma$,
denoted by $N(\mu_{x},\sigma_{x})$. Maximizing
$\mu_{x}-\frac{1}{2}a\sigma_{x}^{2}$ leads to the optimization problem:

$w^{*}_{N}=\arg\max_{w}[w(\mu-r_{0})-\frac{a}{2}w^{2}\sigma^{2}]$ (7)

and the optimal weight $w_{N}$ is given by the Markowitz formula:

$w^{*}_{N}=\frac{\mu-r_{0}}{a\sigma^{2}}$ (8)

Note that weight $w^{*}_{N}$ grows linearly for large $\mu$ without
saturation.

###  4.2 Univariate ALD returns

To compute the expected utility function with ALD, we substitute its PDF into
Eq. 2. The resulting explicit formula is cumbersome, but the optimization-
relevant part can be extracted using the following formal analogy:

$\displaystyle
w^{*}=\arg\max_{w}E_{P}[U_{a}(x,a)]\equiv\left.\arg\max_{w}\big{(}-E_{P}[e^{itx}]\right|_{t=ia}\big{)}$
(9)

The last term is a characteristic function (CF) $E_{P}[e^{itx}]$ of the
probability distribution $P$, well known for most statistical distributions.
For the ALD the CF is given by [7]:

$\phi(t;\mu,\sigma,\mu_{a})=\frac{\exp\left(i\mu
t\right)}{1+\frac{\sigma^{2}t}{2}-i\mu_{a}t},\,\,\,\mu_{a}=\frac{\sigma}{\sqrt{2}}(\frac{1}{\kappa}-\kappa)$
(10)

and results in the following mean $E[X]=\mu+\mu_{a}$ and variance
$Var[X]=\sigma^{2}+\mu_{a}^{2}$. For any random variable $X\sim
ALD(\mu,\sigma,\kappa)$ the following property holds $c+X\sim
ALD(c+\mu,\sigma,\kappa)$. Thus, assuming the outcome of investment is given
by Eq. 6, the portfolio return is distributed according to $ALD(1+r_{0}+w(\mu-
r_{0}),w\sigma,\kappa)$. The optimal weight $w^{*}_{ALD}$ is given by:

$w^{*}_{ALD}=\arg\max_{w}\frac{-e^{-aw(\mu-
r_{0})}}{1-\frac{a^{2}w^{2}\sigma^{2}}{2}+aw\mu_{a}},\,\,\,\,\,s.t.\,\,\frac{a^{2}w^{2}\sigma^{2}}{2}-aw\mu_{a}<1$
(11)

Taking the logarithm, Eq. 11 can be rewritten as follows:

$w^{*}_{ALD}=\arg\max_{w}\left[w(\mu-
r_{0})+\frac{1}{a}\log\left(1-\frac{a^{2}}{2}w^{2}\sigma^{2}+aw\mu_{a}\right)\right],\,s.t.\,\,\frac{a^{2}w^{2}\sigma^{2}}{2}-aw\mu_{a}<1$
(12)

In the symmetrical case $\kappa=1$, the constraint condition simplifies to
$|w|<\frac{\sqrt{2}}{\sigma a}$. If $a\rightarrow 0$, we arrive to the mean-
variance objective function by Taylor expansion.

In order to avoid trivial cluttering, we assume later in the text that an
investor chooses between stock and cash ($r_{0}$=0). Minimization of the
objective function in Eq. 12 leads to the quadratic equation on weights:

$a^{2}\kappa\mu\sigma^{2}w^{2}+a\sigma\left(\sqrt{2}\left(\kappa^{2}-1\right)\mu+2\kappa\sigma\right)w+\sqrt{2}\left(\kappa^{2}-1\right)\sigma-2\kappa\mu=0$
(13)

which has the following (positive) solution:

$w^{*}_{ALD}=\frac{\sqrt{2\left(\kappa^{2}+1\right)^{2}\mu^{2}+4\kappa^{2}\sigma^{2}}-\sqrt{2}\left(\kappa^{2}-1\right)\mu-2\kappa\sigma}{2a\kappa\mu\sigma}$
(14)

In the limit of $\frac{\mu}{\sigma}\to 0$, Eq. 14 simplifies to

$w^{*}_{ALD}=\frac{1}{a\sigma}\left(-\frac{\kappa^{2}-1}{\sqrt{2}\kappa}+\frac{\left(\kappa^{2}+1\right)^{2}}{4\kappa^{2}}\frac{\mu}{\sigma}\right)+O\left(\frac{\mu}{\sigma}\right)^{2}$
(15)

Additionally, if $\kappa\to 1$, we get

$w^{*}_{ALD}=-\sqrt{2}\frac{(\kappa-1)}{a\sigma}+\frac{\mu}{a\sigma^{2}}+O\left(\frac{\mu}{\sigma}\right)^{2}$
(16)

The second item is the standard mean-variance weight from Eq. 8. In the
presence of long only constraint $w>0$, the solution only allows allocation to
the risky asset if $\mu>\sqrt{2}(\kappa-1)\sigma$. For SPX, the monthly
skewness $\kappa_{M}=1.13$, and we obtain
$\mu_{M}>\sqrt{2}(\kappa_{M}-1)\sigma_{M}=0.18*\sigma_{M}$.

Interestingly, the first term in Eq. 15 associated with the skewness of the
distribution appears to be small, but not negligible, relative to the second
term, which resembles the Markowitz solution from Eq. 8. Taking parameters for
S&P 500 daily returns from Table 1, we get

$-\frac{\kappa^{2}-1}{\sqrt{2}\kappa}=-0.06,\,\,\,\frac{\left(\kappa^{2}+1\right)^{2}}{4\kappa^{2}}\frac{\mu}{\sigma}=0.17$
(17)

We see that even a very small skewness $\kappa=1.042$ adds a contribution
comparable to the mean-variance contribution.

The large $\frac{\mu}{\sigma}$ expansion leads to the following asymptotic
solution for optimal portfolio weights:

$w^{*}_{ALD}=\frac{\sqrt{2}}{a\sigma\kappa}+O\left(\frac{1}{\mu}\right)$ (18)

As can be seen, the weights in Eq. 18 do not grow indefinitely, as in the
Markowitz solution 8, but saturate and are independent of the strength of the
signal. This allocation rule can be relevant for quantitative strategies with
a high signal-to-noise ratio.

###  4.3 Marginalization over location parameter $\mu$

The location parameter $\mu$ is often not predictable. To model it, we assume
that $\mu$ follows a normal distribution $\mu\sim N(\mu_{0},\sigma_{0})$ with
parameters $\mu_{0}$ and $\sigma_{0}$, and marginalize (integrate) Eq. 2 over
location parameter $\mu$. The optimal weight $w^{*}$ is given by:

$w^{*}=\arg\max_{w}\int dx\,\,\int
d\mu\,\,U(x,a)P(x;\mu,\sigma)N(\mu;\mu_{0}.\sigma_{0})$ (19)

where $P(x)$ is the distribution of portfolio returns $x$.

To derive the optimal weight for both Gaussian and ALD returns, we rely on the
following identity:
$\int_{-\infty}^{\infty}d\mu\,e^{-aw\mu}\,N(\mu;\mu_{0},\sigma_{0})=e^{\frac{1}{2}a^{2}\sigma_{0}^{2}w^{2}-a\mu_{0}w}$.

For normally distributed returns, the optimal weight $w^{*}_{N}$ after
marginalization over $\mu$ is given by:

$w^{*}_{N}=\arg\max_{w}\left[(-1)e^{-a(\mu_{0}-r)w+\frac{a^{2}}{2}w^{2}(\sigma^{2}+\sigma_{0}^{2})}\right]$
(20)

Taking logarithm of Eq. 20, we arrive to:

$w^{*}_{N}=\arg\max_{w}\left[w(\mu_{0}-r_{0})-\frac{a}{2}w^{2}(\sigma^{2}+\sigma_{0}^{2})\right]=\frac{\mu_{0}-r_{0}}{a(\sigma^{2}+\sigma_{0}^{2})}$
(21)

For ALD returns, the optimal weight $w^{*}$ is:

$w^{*}_{ALD}=\arg\max_{w}\frac{-e^{\frac{1}{2}a^{2}\sigma_{0}^{2}w^{2}-a(\mu_{0}-r)w}}{1-\frac{a^{2}w^{2}\sigma^{2}}{2}+aw\mu_{a}},\,\,\,\,\,s.t.\,\,\frac{a^{2}w^{2}\sigma^{2}}{2}-aw\mu_{a}<1$
(22)

The optimal weight is given by a solution of cubic equation with respect to
$w$:

$w^{3}a^{3}\kappa\sigma^{2}\sigma_{0}^{2}+w^{2}a^{2}\left(\sqrt{2}\left(\kappa^{2}-1\right)\sigma\sigma_{0}^{2}-\kappa\mu_{0}\sigma^{2}\right)-$
(23)
$-wa\left(\sqrt{2}\left(\kappa^{2}-1\right)\mu_{0}\sigma+2\kappa(\sigma^{2}+\sigma_{0}^{2})\right)-\sqrt{2}\left(\kappa^{2}-1\right)\sigma+2\kappa\mu_{0}=0$

The roots of the cubic equation can be found using Cardano’s formula. Of the
three solutions, only one has a negative second derivative with respect to
$w$, which guarantees a maximum required in Eq. 22. Furthermore, we do not
provide a cumbersome formula for the exact solution, but study the three most
important asymptotics. We assume that $\kappa\to 1$ in all cases.

For investors with low confidence in the predicted expected return
$\frac{\mu_{0}}{\sigma_{0}}\to 0$, the optimal weight is given by:

$\lim_{\frac{\mu_{0}}{\sigma_{0}}\to
0}w^{*}=\frac{mean}{a(\sigma^{2}+\sigma_{0}^{2})}+O(\kappa-1)^{2},\,\,mean=\mu_{0}+\sqrt{2}(1-\kappa)\sigma$
(24)

For investors with high confidence in the predicted expected return
$\frac{\mu_{0}}{\sigma_{0}}\to\infty$, the optimal weight is given by:

$\lim_{\frac{\mu_{0}}{\sigma_{0}}\to\infty}w^{*}=\frac{\sqrt{2}}{a\sigma\kappa}+O(\kappa-1)^{2}\approx\frac{\sqrt{2}}{a\sigma}$
(25)

The limit of $\sigma_{0}\to 0$ is given by:

$\lim_{\sigma_{0}\to
0}w^{*}=-\frac{\sqrt{2}(\kappa-1)}{a\sigma}+\frac{\sqrt{2\mu_{0}^{2}+\sigma^{2}}-\sigma}{a\mu_{0}\sigma}+O(\kappa-1)^{2}$
(26)

##  5 Multivariate case

In the multivariate case of $N$ assets, the vectors
$\boldsymbol{w},\boldsymbol{\mu},\boldsymbol{\mu}_{a}\in R^{N}$ and covariance
matrix $\boldsymbol{\Sigma}\in R^{N\times N}$. In the case of multivariate
Gaussian returns $\boldsymbol{r}_{t}\sim
N(\boldsymbol{\mu},\boldsymbol{\Sigma})$, the expected utility function is
proportional to:

$E_{N}[U_{a}]\sim
exp(a\boldsymbol{\mu}^{T}\boldsymbol{w}-\frac{1}{2}a^{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w})$
(27)

After taking the logarithm of $E_{N}[U_{a}]$, the optimization problem
becomes:

$\boldsymbol{w}^{*}=\arg\max_{\boldsymbol{w}}\left[\boldsymbol{\mu}^{T}\boldsymbol{w}-\frac{a}{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}\right]$
(28)

which has the following solution:

$\boldsymbol{w}^{*}=\frac{1}{a}\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}$ (29)

Furthermore, we find allocation rules for portfolio returns distributed
according to the multivariate ALD. The distribution is defined by its location
vector $\boldsymbol{\mu}$, scale matrix $\mathbf{\Sigma}$, and asymmetry
parameter $\boldsymbol{\mu_{a}}$. The CF of multivariate ALD is given by [7]:

$\phi(\boldsymbol{t};\boldsymbol{\mu},\boldsymbol{\Sigma},\boldsymbol{\mu}_{a})=\frac{\exp\left(i\boldsymbol{\mu}^{T}\boldsymbol{t}\right)}{1+\frac{1}{2}\boldsymbol{t}^{T}\boldsymbol{\Sigma}\boldsymbol{t}-i\boldsymbol{\mu}_{a}^{T}\boldsymbol{t}}$
(30)

Here, $\boldsymbol{\mu}_{a}$ is still a scale-dependent quantity. There are
multiple ways to make it scale-independent. We define $\boldsymbol{\mu}_{a}$
using scale-independent parameters $\boldsymbol{\kappa}$ as follows:

$\boldsymbol{\mu}_{a}=\frac{1}{\sqrt{2}}\sqrt{Diag(\boldsymbol{\Sigma})}(\boldsymbol{\kappa}^{-1}-\boldsymbol{\kappa})$
(31)

The mean and variance are given by:

$Mean=\boldsymbol{\mu}+\boldsymbol{\mu}_{a},\,\,Var=\boldsymbol{\Sigma}+\boldsymbol{\mu}_{a}^{T}\boldsymbol{\mu}_{a}$
(32)

The expected utility function $E_{ALD}[U_{a}]$ is proportional to:

$E_{ALD}[U_{a}]\sim\frac{-\exp\left(-a\boldsymbol{\mu}^{T}\boldsymbol{w}\right)}{1-\frac{a^{2}}{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}+a\boldsymbol{\mu}_{a}^{T}\boldsymbol{w}}$
(33)

Taking the logarithm of $E_{ALD}[U_{a}]$ and omitting irrelevant constants, we
have the following optimization problem:

$\boldsymbol{w}^{*}=\arg\max_{\boldsymbol{w}}\left[\boldsymbol{\mu}^{T}\boldsymbol{w}+\frac{1}{a}\log\left(1-\frac{a^{2}}{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}+a\boldsymbol{\mu}_{a}^{T}\boldsymbol{w}\right)\right]$
(34)

such that
$\frac{a^{2}}{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}-a\boldsymbol{\mu}_{a}^{T}\boldsymbol{w}<1$.
Maximization of Eq. 34 results in the following system of quadratic equations
with respect to $\boldsymbol{w}$:

$\boldsymbol{\mu}\left(1-\frac{a^{2}}{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}+a\boldsymbol{\mu}_{a}^{T}\boldsymbol{w}\right)+(\boldsymbol{\mu}_{a}-a\boldsymbol{\Sigma}\boldsymbol{w})=0$
(35)

For the symmetric case with $\mu_{a}=0$, the solution is given by:

$\boldsymbol{w}^{*}=\frac{1-\frac{a^{2}}{2}d}{a}\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}$
(36)

where $d=\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}$. Substituting
the solution (36) into the definition of $d$, we obtain a quadratic equation
with the following solution:

$d=\frac{2(1+q+\sqrt{1+2q})}{a^{2}q}$ (37)

where $q=\boldsymbol{\mu}^{T}\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}$. In the
general case, Eq. 35 can be solved numerically. The solve_poly_system solver
in the Sympy Python package is a reasonable choice for this task.

Similar to the univariate case considered in Section 4.3, we study the
expected utility function marginalized over $\boldsymbol{\mu}$ with
$\boldsymbol{\mu}\sim N(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{0})$. Here,
$\boldsymbol{\Sigma}_{0}$ is a diagonal matrix whose elements are equal to the
variance of individual components of $\boldsymbol{\mu}$. To derive the optimal
weight for both multivariate Gaussian and ALD returns, we rely on the
following identity:
$\int_{-\infty}^{\infty}d\boldsymbol{\mu}\,e^{-a\boldsymbol{w}^{T}\boldsymbol{\mu}}\,N(\boldsymbol{\mu};\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{0})=e^{\frac{1}{2}a^{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}_{0}\boldsymbol{w}-a\boldsymbol{\mu}_{0}^{T}\boldsymbol{w}}$,

For normally distributed returns, the optimal weights $\boldsymbol{w}^{*}$
after marginalization over $\boldsymbol{\mu}$ are given by:

$\boldsymbol{w}^{*}=\arg\max_{\boldsymbol{w}}\left[-e^{\frac{1}{2}a^{2}\boldsymbol{w}^{T}\boldsymbol{(}\boldsymbol{\Sigma}+\boldsymbol{\Sigma}_{0})\boldsymbol{w}-a\boldsymbol{\mu}_{0}^{T}\boldsymbol{w}}\right]$
(38)

This equation can be rewritten as follows:

$\boldsymbol{w}^{*}=\arg\max_{\boldsymbol{w}}\left[\boldsymbol{\mu}_{0}^{T}\boldsymbol{w}-\frac{a}{2}\boldsymbol{w}^{T}\boldsymbol{(}\boldsymbol{\Sigma}+\boldsymbol{\Sigma}_{0})\boldsymbol{w}\right]$
(39)

The optimal weights are given by:

$\boldsymbol{w}^{*}=\frac{1}{a}(\boldsymbol{\Sigma}+\boldsymbol{\Sigma}_{0})^{-1}\boldsymbol{\mu}_{0}$
(40)

The optimal weights $\boldsymbol{w}^{*}$ for ALD distributed returns are given
by:

$\boldsymbol{w}^{*}=\arg\max_{\boldsymbol{w}}\frac{-e^{\frac{1}{2}a^{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}_{0}\boldsymbol{w}-a\boldsymbol{\mu}_{0}^{T}\boldsymbol{w}}}{1-\frac{a^{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}}{2}+a\boldsymbol{\mu}_{a}^{T}\boldsymbol{w}},\,\,s.t.\,\,\frac{a^{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}}{2}-a\boldsymbol{\mu}_{a}^{T}\boldsymbol{w}<1$
(41)

This equation can be rewritten as follows:

$\boldsymbol{w}^{*}=\arg\max_{\boldsymbol{w}}\left[\boldsymbol{\mu}_{0}^{T}\boldsymbol{w}-\frac{a}{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}_{0}\boldsymbol{w}+\frac{1}{a}\log\left(1-\frac{a^{2}}{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}+a\boldsymbol{\mu}_{a}^{T}\boldsymbol{w}\right)\right]$
(42)

such that
$\frac{a^{2}}{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}-a\boldsymbol{\mu}_{a}^{T}\boldsymbol{w}<1$.
Eq. 42 can be solved by a convex numerical optimizer such as CVXPY [9] even
with additional constraints such as long-only, maximum number of position,
turnover, and other constraints.

Unconstrained solution of Eq. 42 can be obtained by taking a derivative with
respect of $\boldsymbol{w}$. Thus, we arrive at the following system of cubic
equations:

$\boldsymbol{(}\boldsymbol{\mu}_{0}-a\boldsymbol{\Sigma}_{0}\boldsymbol{w})\left(1-\frac{a^{2}}{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}+a\boldsymbol{\mu}_{a}^{T}\boldsymbol{w}\right)+(\boldsymbol{\mu}_{a}-a\boldsymbol{\Sigma}\boldsymbol{w})=0$
(43)

For the symmetric case with $\boldsymbol{\mu}_{a}=0$, the solution is given
by:

$\boldsymbol{w}^{*}=\frac{(1-\frac{a^{2}}{2}d)\boldsymbol{\mu}_{0}}{a(\boldsymbol{\Sigma}+(1-\frac{a^{2}}{2}d)\boldsymbol{\Sigma}_{0})}$
(44)

where $d=\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}$. Substituting
the solution (44) into the definition of $d$, we obtain a nonlinear equation
for d which can be solved numerically. Note that the role of uncertainty
$\boldsymbol{\Sigma}_{0}$ of the location parameter $\boldsymbol{\mu}$ is
similar to the regularization (shrinkage) approach [10], but has a clearer
interpretation. Moreover, there is no need to estimate a shrinkage constant.
In the general case, Eq. 43 can be solved numerically.

##  6 Worst-case optimization

One of the most prominent trends in trading is the increased use of leverage,
both in stocks and in derivatives such as futures and options. Leverage
enables traders to control a large contract value with a relatively small
amount of capital. Leveraged trading often carries a non-negligible
probability of ruin, as it amplifies both potential gains and losses.
Consequently, to manage such risks, asset allocation under the worst-case
outcome scenario often emerges as the optimal strategy in trading. In worst-
case optimization, we maximize the expected value of the utility function
under the worst-case outcome for the expected returns:

$\boldsymbol{w}^{*}=\arg\max_{\boldsymbol{w}}\min_{\boldsymbol{\mu}\in D}E[U]$
(45)

where $D$ is the worst outcome under a model $\boldsymbol{\mu}\sim
N(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{0})$.

###  6.1 Worst-case optimization for the risk-neutral case

If the risk component is neglected, the optimization problem is given by:

$\boldsymbol{w}^{*}=\arg\max_{\boldsymbol{w}}\min_{D}\boldsymbol{\mu}^{T}\boldsymbol{w},\,\,s.t.\sum_{i=1}^{N}w_{i}=1,\,w_{i}\geq
0$ (46)

We assume that the expected returns follow a multivariate normal distribution
with mean $\boldsymbol{\mu}_{0}$ and diagonal variance
$\boldsymbol{\Sigma}_{0}$. Given that $\boldsymbol{\mu}\sim
N(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{0})=\boldsymbol{\mu}_{0}+\boldsymbol{Y}$
and $\boldsymbol{Y}\sim N(0,\boldsymbol{\Sigma}_{0})$, the worst outcome for
the weighted sum is given by:

$\min_{D}\boldsymbol{\mu}^{T}\boldsymbol{w}=\boldsymbol{\mu}_{0}^{T}\boldsymbol{w}+\min\boldsymbol{Y}*\max_{i}w_{i}$
(47)

where $\max_{i}w$ is the largest weight from $\boldsymbol{w}$. We estimate the
minimal value of error $\min\boldsymbol{Y}$ under the assumption that
$\boldsymbol{Y}\sim N(0,\boldsymbol{\Sigma}_{0})$. Because matrix
$\boldsymbol{\Sigma}_{0}$ is diagonal, all $Y_{i}$ are independent, and the
problem is to find the expected value $\langle Y_{min}\rangle$ of the minimum
of N random variables $Y_{i}\sim N(0,\sigma^{0}_{i})$ with Gaussian CDF
$F_{i}(Y_{i}\leq t)$. The CDF of the minimum of such variables is given by:

$F_{Y}(y)=P(\min_{i}Y_{i}\leq t)=1-\prod_{i}^{N}(1-F_{i}(Y_{i}\leq t))$ (48)

The PDF and the expected value of the minimum are given by:

$f_{Y}(y)=\frac{dF_{Y}(y)}{dy}=\sum_{i=1}^{N}f_{i}(y)*\prod_{i\neq
i}(1-F_{j}(y)),\,\,\langle
Y_{min}\rangle=\int_{-\infty}^{-\infty}dy\,yf_{Y}(y)$ (49)

the expected value $\langle Y_{min}\rangle$ can be evaluated numerically by
substituting into Eq. 49 the probability density function $f(x)$ and
cumulative distribution function $F(x)$ of a normal distribution with mean
$\mu$ and variance $\sigma^{2}$:

$f(x)=\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}},\,\,\,F(x)=\frac{1}{2}\left[1+\operatorname{erf}\left(\frac{x-\mu}{\sigma\sqrt{2}}\right)\right]$
(50)

where $\operatorname{erf}(x)$ is the error function. In the case of equal
$\sigma^{0}_{i}=\tilde{\sigma}_{0}$ and large N, the $\langle Y_{min}\rangle$
has the form $\langle Y_{min}\rangle\approx-\tilde{\sigma}_{0}\sqrt{2\log N}$.

Eq. 47 is a convex problem, and any solution satisfying the Karush-Kuhn-Tucker
(KKT) conditions can be expected to be optimal. The KKT method generalizes the
Lagrange multiplier method in the presence of inequality. When the objective
function and constraints are convex, the KKT conditions are not only necessary
but also sufficient for global optimality.

As we show in Appendix A, the solution of Eq. 47, derived with Lagrange
multipliers and KKT conditions, exhibits a simple behavior. For small $\langle
Y_{min}\rangle$, all weights are allocated to an asset with the largest
$\mu^{0}_{i}$. As $\langle Y_{min}\rangle$ starts to increase, the allocations
change to a uniform allocation for $k$ assets with the largest $\mu^{0}_{i}$,
where $k$ is an implicit function of $\boldsymbol{\mu}^{0}$ and
$\boldsymbol{\sigma}^{0}$. For large $\langle Y_{min}\rangle$, under the
condition $|\langle
Y_{min}\rangle|/N>\left(\langle\boldsymbol{\mu}^{0}\rangle-\min_{i}\mu^{0}_{i}\right)$,
the solution is an equal weights $\frac{1}{N}$ allocation. This condition is
satisfied early if the difference between the average
$\langle\boldsymbol{\mu}^{0}\rangle$ and the minimal expected return
$\min_{i}\mu^{0}_{i}$ is small, which typically occurs for similar assets. The
solution can be understood qualitatively, as the loss function in Eq. 47 is
dominated by the $\max_{i}w$ term, which penalizes any allocation deviating
from a uniform $1/N$. The $\max_{i}w$ term can also be interpreted as an
infinity norm of the vector $\boldsymbol{w}$, given the following identity:
$\max_{i}w=|w|_{\infty}$.

###  6.2 Worst-case optimization with risk for the Gaussian return model

If returns are distributed according to a multivariate normal distribution, we
optimize Eq. 28 with the worst expected returns from Eq. 47:

$\boldsymbol{w}^{*}=\arg\max_{w}\left[\boldsymbol{\mu}_{0}^{T}\boldsymbol{w}+\min\boldsymbol{Y}*\max_{i}w_{i}-\frac{a}{2}\boldsymbol{w}^{T}\Sigma\boldsymbol{w}\right],\,\,s.t.\sum_{i=1}^{N}w_{i}=1,\,w_{i}\geq
0$ (51)

The solution can be obtained KKT conditions. For most practical cases,
numerical optimization with CVXPY can provide a globally optimal solution with
minimal effort.

Eq. 51 has an important limiting case when the contribution
$\boldsymbol{\mu}_{0}^{T}\boldsymbol{w}$ can be neglected. In this case, we
have:

$\boldsymbol{w}^{*}=\arg\min_{w}\left[\max_{i}w_{i}+\frac{b}{2}\boldsymbol{w}^{T}\Sigma\boldsymbol{w}\right],\,\,s.t.\sum_{1}^{N}w_{i}=1,\,w_{i}\geq
0,\,\,b>0$ (52)

where $b=\frac{a}{|\min\boldsymbol{Y}|}$. Because both the terms are convex,
the optimization problem is convex, which ensures the existence of a unique
global optimum. We discuss KTT solution of Eq. 52 in Appendix B.

It is instructive to analyze two terms in the objective function Eq. 52
separately. The first term, $\max_{i}w_{i}$, penalizes the maximum weight in
the portfolio. Minimizing this term alone would lead to equal-weighted or more
evenly distributed portfolio allocations. The second term,
$\boldsymbol{w}^{T}\Sigma\boldsymbol{w}$, can be thought of as a measure of
portfolio variance. Minimizing this term alone would lead to a minimum
variance portfolio. When both terms are combined, the optimization problem
aims to balance the trade-off between minimizing portfolio variance and
avoiding concentrated portfolio allocations. As the constant term increases,
the optimization problem will put more emphasis on avoiding concentrated
allocations, leading to a more evenly distributed portfolio.

A portfolio with similar characteristics is a risk parity(RP) portfolio, which
is also known to provide a trade-off between the minimum variance and equally
weighted portfolios [11]. In the case of equal correlation, RP allocation is
equal volatility allocation. The RP objective function is not convex, and to
compute the optimal weights, nonlinear optimization methods, such as
sequential quadratic programming, are required. Most importantly, the RP
approach misses the critical ingredient: any information about expected
returns and their associated uncertainties.

Our formulation offers a convex alternative to RP and includes an adjustable
parameter $b$ that governs the trade-off between portfolio volatility risk and
the worst-case scenario for the expected return prediction. The information
about expected return mean $\boldsymbol{\mu}_{0}$ can also be included. We
refer to this class of portfolios as minimax (MM) portfolios.

We propose a heuristic rule to find a starting value for parameter $b$ that
balances the lack of excessive concentration and correlation structure of a
portfolio. Let $S$ be the entropy of portfolio weights
$S=-\sum_{i=1}^{N}w_{i}\log w_{i}$ and $N_{eff}$ is the effective number of
securities in a portfolio $N_{eff}=\exp{S}$. A plausible $b^{*}$ is given by
the condition that the entropy of the MM portfolio is equal to the average of
entropy for the equal weight portfolio $S_{eq}$ and the minimum variance
portfolio $S_{mv}$:

$S_{b^{*}}=\frac{1}{2}(S_{eq}+S_{mv})$ (53)

This rule provides a starting point for the parameter $b$, which can be
further refined through numerical simulations or other optimization
techniques. Alternatively, an investor can choose $b^{*}$ to obtain a specific
value of $N_{eff}(b^{*})/N$. $N_{eff}=N$ for an equally weighted portfolio and
$N_{eff}=1$ for a single stock portfolio.

###  6.3 Worst case optimization with risk for the ALD return model

Assuming that returns are distributed according to a multivariate ALD, we
optimize Eq. 34 with the worst expected returns from Eq. 47:

$\boldsymbol{w}^{*}=\arg\max_{\boldsymbol{w}}\left[\boldsymbol{\mu}_{0}^{T}\boldsymbol{w}+\min
Y*\max_{i}w_{i}+\frac{1}{a}\log\left(1-\frac{a^{2}}{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}+a\boldsymbol{\mu}_{a}^{T}\boldsymbol{w}\right)\right]$
(54) $s.t.\sum_{i=1}^{N}w_{i}=1,\,w_{i}\geq
0,\,\,\,\,\frac{a^{2}}{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}-a\boldsymbol{\mu}_{a}^{T}\boldsymbol{w}<1$

For general combinations of parameters, numerical optimization with CVXPY can
provide a global optimum for Eq. 54 in minimal time.

If the contribution of expected return
$\boldsymbol{\mu}_{0}^{T}\boldsymbol{w}$ is small, the optimization function
can be rewritten as follows:

$\boldsymbol{w}^{*}=\arg\min_{w}\left[\max_{i}w_{i}-b\log\left(1-\frac{a^{2}}{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}+a\boldsymbol{\mu}_{a}^{T}\boldsymbol{w}\right)\right]$
(55) $s.t.\sum_{i=1}^{N}w_{i}=1,\,w_{i}\geq
0,\,\,\,\,\frac{a^{2}}{2}\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}-a\boldsymbol{w}^{T}\boldsymbol{\mu}_{a}<1$

Empirically, we observe that the solution of Eq. 55 is often close to the one
defined by the Gaussian case in Eq. 52. However, for large $N$, we needed to
switch from the ECOS optimizer to a more sophisticated SCS in CVXPY.

##  7 Imposing block structure on a covariance matrix

Covariance and correlation coefficients are frequently used to measure the
relationship between stocks. A high correlation coefficient may arise either
because the stocks are directly correlated with one another or because they
are correlated with a similar index through their beta values. In the latter
case, the index acts as a confounder, and when controlling for index returns,
the stocks may move independently. Mathematically, this means that the stocks
are conditionally independent, given the return of the index. Covariance is
subject to numerical noise arising from finite samples and the existence of
confounding variables. Moreover, spurious correlations are abundant in
financial time series [12].

Although machine learning algorithms, such as graphical lasso [13], can be
used to learn the true conditional dependence (graphical) structure,
experienced portfolio managers likely have this structure in mind before
allocating real capital to an investment strategy. As a first approximation,
the information on conditional independence can be incorporated into portfolio
construction by utilizing a block diagonal structure derived from shared
characteristics, such as sector, investment style, geographical region, or
asset class. This structure assumes a stronger relationship within blocks and
weaker or no causal relationship between blocks, simplifying portfolio risk
estimation and allocation. As a result, many spurious correlations between
blocks disappear. Our objective is to estimate the covariance and precision
matrices with a superimposed structure (adjacency) matrix $\boldsymbol{A}$.

In this section, we examine how to impose a graphical structure on covariance
or precision matrices, which often results in a block diagonal structure.
Although a component of the covariance matrix $\boldsymbol{\Sigma}$ represents
marginal dependence between corresponding random variables, elements of the
precision matrix $\boldsymbol{\Theta}$ signify their conditional independence,
providing information about partial correlations between each pair of
variables after accounting for the effects of all other variables. The
elements of the partial correlation matrix $\varrho_{ij}$ can be obtained from
the precision matrix $\Theta$ as follows:
$\varrho_{ij}=-\frac{\Theta_{ij}}{\sqrt{\Theta_{ii}\Theta_{jj}}}$. The off-
diagonal elements of $\varrho$ can be used as an adjacency matrix of the
portfolio graphical structure.

Consider a portfolio consisting of two energy sector stocks, $E_{1}$ and
$E_{2}$, and two healthcare stocks, $H_{1}$ and $H_{2}$. The returns of stocks
in different sectors are conditionally independent; if we are interested in
the return of $E_{1}$ and know the returns of $E_{2}$, the information on
healthcare stocks returns becomes redundant. In other words, the block
diagonal structure implies that there is no linear relationship between the
assets when controlling for the other assets in the portfolio.

We demonstrate the concept using the Graphical Sparse Precision Matrix
Estimation (GSPME) method, as implemented in the GraphSPME package [14]. This
approach combines the precision matrix estimation routine with respect to a
graph [15] with the automatic and adaptive shrinkage of sample covariance
matrices developed in [16]. Importantly, the resulting precision matrix
estimates are guaranteed to be positive semi-definite, which allows for the
utilization of efficient matrix factorization methods, such as Cholesky
decomposition, in order to perform matrix inversion and subsequently obtain a
covariance matrix.

As an example, we consider a portfolio consisting of Exxon Mobil (XOM) and
Chevron (CVX) energy stocks, and Pfizer (PFE) and Merck (MRK) healthcare
stocks, using daily returns for the year 2022. Assuming conditional
independence between the energy and healthcare sectors, the structure
information about conditional independence is encoded by an adjacency matrix
$A$, which contains zeros and ones:

$\boldsymbol{A}=\left(\begin{tabular}[]{rrrr}1&1&0&0\\\ 1&1&0&0\\\ 0&0&1&1\\\
0&0&1&1\\\ \end{tabular}\right)$

We compare covariance $\boldsymbol{\Sigma}$ and precision matrices
$\boldsymbol{\Theta}$ and correlation $\boldsymbol{C}orr$ and partial
correlation $\boldsymbol{\varrho}$ matrices using direct calculation and the
proposed method for the portfolio above in Appendix C.

The problem of inverting a covariance matrix becomes severe when the number of
assets $N$ is larger than the number of observations $d$ as the matrix becomes
singular. We explored this problem using S&P 500 stocks, with an adjacency
matrix $\boldsymbol{A}$ given by sector constraints, using the daily returns
of continuously traded stocks between 2010 and 2022. For ten industry sectors
(excluding real estate), we have twelve yearly samples with $d=252$ trading
days per year (on average) and a total of $N=354$ stocks. The model precision
matrix $\boldsymbol{\Theta}_{model}$ for each year is full rank, has a block
diagonal structure with the largest block size much smaller than the number of
observations in a year ($d=252$), is semi-positive definite by construction,
and allows for fast inversion and power operations. The condition number of
the estimated covariance matrix $\boldsymbol{\Theta}_{model}^{-1}$ is only
2,134 for the year 2022. In summary, we believe that the precision matrix
$\boldsymbol{\Theta}$ is a primary object for problems that require block
(top-down) constraints, and the covariance matrix $\boldsymbol{\Sigma}$ can be
derived from $\boldsymbol{\Theta}$ through inversion operations.

Another source of noise in covariance matrix estimation is the non-
stationarity of return time series. The noise encompasses two components:
variance and correlation structure. In order to mitigate these effects, it is
prudent to employ a separation strategy as proposed in [17], where standard
deviations and correlations are modeled independently. The covariance is
decomposed as follows:

$\boldsymbol{\Sigma}=\boldsymbol{\Lambda}\boldsymbol{R}\boldsymbol{\Lambda}$
(56)

where, $\boldsymbol{\Lambda}$ represents a diagonal matrix containing
volatility elements $\sigma_{i}$ and $R$ denotes the corresponding correlation
matrix. For volatility estimation in $\boldsymbol{\Lambda}$, one can utilize
forward-looking volatility estimations derived from option markets. The
correlation matrix $\boldsymbol{R}$ remains relatively stable, except during
regime changes such as market crashes when all assets become highly correlated
and move in unison.

##  8 Portfolio allocation for long-term investors

The long-term cross-sectional returns of major stock indexes are well
approximated by log-normal distribution [5]. The expected utility function is
given by substitution of the CRRA utility $U_{r}(x,\gamma)$ and log-normal
distribution function $P(x)=LN(x;\mu,\sigma)$ into Eq. 2:

$\displaystyle\begin{split}&E[U_{r}(x,\gamma)]=\int
dx\,\,U_{r}(x,\gamma)\,LN(x;\mu,\sigma)=\frac{e^{(1-\gamma)\mu+(1-\gamma)^{2}\frac{\sigma^{2}}{2}}-1}{1-\gamma}\\\
&E[U_{r}(x,\gamma=1)]=\mu\end{split}$ (57)

The maximization of $E[U_{r}(x,\gamma)]$ has a clear interpretation for
$\gamma=0$ and $\gamma=1$. In the former case, the utility function is linear
(risk-neutral) $U_{r}(x,0)=x-1$ and an investor should maximize the average
return of the portfolio:

$w^{*}_{LN}=\max e^{\mu+\frac{\sigma^{2}}{2}}=\max E[x]$ (58)

In the latter case, the utility function is logarithmic $U_{r}(x,1)=\ln(x)$,
and an investor should maximize the median return of the portfolio:

$w^{*}_{LN}=\max e^{\mu}=\max median[x]$ (59)

Assume that each stock price $X^{i}$ in a portfolio follows a geometric
Brownian motion:

$\frac{dX^{i}_{t+1}}{X^{i}_{t}}=\mu_{i}\,dt+\sigma\,dW_{t},$ (60)

with the percentage drift distributed according to the normal distribution
$\mu_{i}\sim N(\mu_{d},\sigma_{d})$. In this model [5], the total return
$\rho_{i}=X^{i}_{T}/X^{i}_{0}$ of a randomly chosen stock $i$ at time $T$
follows a log-normal distribution with
$\mu_{m}=\mu_{d}T-\frac{1}{2}\sigma^{2}T$ and
$\sigma_{m}^{2}=\sigma^{2}T+\sigma_{d}^{2}T^{2}$. The mean and median of the
return distribution are given by:

$mean[\rho]=e^{\mu_{m}+\frac{\sigma_{m}^{2}}{2}}=e^{\mu_{d}T+\frac{1}{2}\sigma_{d}^{2}T^{2}}$
(61) $median[\rho]=e^{\mu_{m}}=e^{\mu_{d}T-\frac{1}{2}\sigma^{2}T}$ (62)

A risk-neutral investor with utility function $U_{r}(x,0)=x-1$ should maximize
the expected value of the drift $\mu_{d}$ and, most importantly, maximize the
dispersion of the drift $\sigma_{d}$ by taking diversified risky bets and
disregarding the volatility of the stock in the portfolio. Drift dispersion
${\sigma_{d}}^{2}T^{2}$ represents the effect of continuous compounding of
winners pushing the average return up, while a large body of distribution is
concentrated around the mode and goes to zero. In turn, an investor with
logarithmic utility $U_{r}(x,1)=\ln(x)$ should maximize the expected value of
the drift $\mu_{d}$ while minimizing the volatility exposure $\sigma^{2}T$.

##  9 Conclusions

This paper introduces a series of models designed to enhance the quality of
asset allocation rules. Analyzing historical stock index returns across
various time scales, we find that the ALD serves as a suitable return model
for daily, weekly, and monthly returns. We derive allocation rules for ALD
returns and systematically compare them to the mean-variance approach. Our
model enables to study the average and the worst-case scenarios, taking into
account the uncertainty associated with the expected returns. We investigate
the integration of block structures into portfolio allocation rules. By
incorporating block constraints into precision matrix estimation, we promote
conditional independence between inter-block constituents and facilitate the
inversion of singular covariance or precision matrices under specific
conditions. Finally, we develop allocation rules for log-normal distributed
cross-sectional returns, laying the groundwork for long-term portfolio
allocations.

## References

  * Markowitz [1952] Harry Markowitz.  Portfolio selection.  _The Journal of Finance_ , 7(1):77–91, 1952.
  * Demiguel et al. [2009] Victor Demiguel, Lorenzo Garlappi, and Raman Uppal.  Optimal versus naive diversification: How inefficient is the 1/n portfolio strategy?  _Review of Financial Studies_ , 22, 05 2009.
  * Taleb [2008] Nassim N. Taleb.  _The Black Swan: The Impact of the Highly Improbable_.  Random House, London, 1 edition, 2008.
  * Cochrane [2021] John H. Cochrane.  Portfolios for long-term investors.  _Review of Finance_ , 26(1):1–42, 2021.
  * Markov and Markov [2023] Maxime Markov and Vladimir Markov.  The impact of big winners on passive and active equity investment strategies, 2022.  URL https://doi.org/10.48550/arXiv.2210.09302.
  * Campbell and Viceira [2001] John Y. Campbell and Luis M. Viceira.  _Strategic asset allocation: Portfolio choice for long-term investors_.  Oxford University Press, 2001.
  * Kotz et al. [2001] Samuel Kotz, Tomasz J. Kozubowski, and Krzysztof Podgorski.  _The Laplace Distribution and Generalizations: A Revisit with New Applications_.  Springer Science and Business Media, 2001.
  * Salvatier et al. [2016] John Salvatier, Thomas V. Wiecki, and Christopher Fonnesbeck.  Probabilistic programming in python using PyMC3.  _PeerJ Computer Science_ , 2:e55, 2016.
  * Diamond and Boyd [2016] Steven Diamond and Stephen Boyd.  CVXPY: A Python-embedded modeling language for convex optimization.  _Journal of Machine Learning Research_ , 17(83):1–5, 2016.
  * Ledoit and Wolf [2003] Olivier Ledoit and Michael Wolf.  Improved estimation of the covariance matrix of stock returns with an application to portfolio selection.  _Journal of Empirical Finance_ , 10(5):603–621, 2003.
  * Maillard et al. [2010] Sébastien Maillard, Thierry Roncalli, and Jerome Teiletche.  The properties of equally weighted risk contribution portfolios.  _The Journal of Portfolio Management_ , 36:60–70, 07 2010\.
  * Leinweber [2007] David Leinweber.  Stupid data miner tricks: Overfitting the $S\&P$ 500\.  _The Journal of Investing_ , 16:15–22, 02 2007.
  * Friedman et al. [2008] Jerome Friedman, Trevor Hastie, and Robert Tibshirani.  Sparse inverse covariance estimation with the graphical lasso.  _Biostatistics (Oxford, England)_ , 9:432–41, 08 2008.
  * Lunde et al. [2022] Berent Anund Stromnes Lunde, Feda Curic, and Sondre Sortland.  Graphspme: Markov precision matrix estimation and asymptotic stein-type shrinkage.  2022\.  URL https://doi.org/10.48550/arXiv.2205.07584.
  * Le and Zhong [2021] Thien Le and Ping‐Shou Zhong.  High‐dimensional precision matrix estimation with a known graphical structure.  _Stat_ , 11, 09 2021.
  * Touloumis [2015] Anestis Touloumis.  Nonparametric stein-type shrinkage covariance matrix estimators in high-dimensional settings.  _Computational Statistics & Data Analysis_, 83:251–261, 03 2015.
  * Barnard et al. [2000] John Barnard, Robert McCulloch, and Xiao Li Meng.  Modeling covariance matrices in terms of standard deviations and correlations, with application to shrinkage.  _Statistica Sinica_ , 10:1281–1311, 2000.

### Appendix A: Worst-case optimization in risk-neutral case

The objective of the worst-case optimization in the risk-neutral case is given
by the following equation:

$\boldsymbol{w}^{*}=\arg\max_{w}\left[\boldsymbol{\mu}_{0}^{T}\boldsymbol{w}-c\cdot
r\right],\,s.t.\,\,w_{i}\geq 0,\sum_{i=1}^{N}w_{i}=1,\,\,w_{i}\leq r$ (63)

We can define the Lagrangian function as follows:

$L(\boldsymbol{w},\lambda,\boldsymbol{\eta},\boldsymbol{\nu})=\boldsymbol{\mu}_{0}^{T}\boldsymbol{w}-c\cdot
r+\lambda\cdot(1-\sum_{i=1}^{N}w_{i})+\sum_{i=1}^{N}\eta_{i}\cdot
w_{i}-\sum_{i=1}^{N}\nu_{i}\cdot(w_{i}-r)$ (64)

KKT conditions:

  1. 1.

Stationarity:
$\nabla_{\boldsymbol{w}}L=\boldsymbol{\mu}_{0}-\lambda\cdot\mathbf{1}+\boldsymbol{\eta}-\boldsymbol{\nu}=0,\,\,\,\,\partial_{r}L=-c+\mathbf{1}^{T}\boldsymbol{\nu}=0$

  2. 2.

Primal feasibility: $w_{i}\geq 0,\sum_{i=1}^{N}w_{i}=1,w_{i}\leq r$

  3. 3.

Dual feasibility: $\eta_{i}\geq 0,\nu_{i}\geq 0$

  4. 4.

Complementary slackness: $\eta_{i}\cdot w_{i}=0,\,\,\nu_{i}\cdot(w_{i}-r)=0$

Let us assume that $\mu_{0}^{i}$ are positive and ordered in descending order.
For non-zero $w_{i}$, we have the following equation:

$\eta_{i}=0,\nu_{i}=\mu_{0}^{i}-\lambda,\nu_{i}\geq 0$ (65)

Using the second stationarity condition, the optimal value of $\lambda^{*}$ is
given by $\sum_{i=1}^{N}(\mu_{0}^{i}-\lambda^{*})^{+}=c$, where
$(x)^{+}=\max(0,x)$. This leads to the following condition for $\lambda^{*}$:
$\mu_{0}^{k+1}<\lambda^{*}<\mu_{0}^{k}$. The number of non-zero weights $k$ is
defined as: $\sum_{i=1}^{k}(\mu_{0}^{i}-\mu_{0}^{k})<c$ and
$\sum_{i=1}^{k+1}(\mu_{0}^{i}-\mu_{0}^{k+1})>c$. Thus,
$\lambda^{*}=\frac{1}{k}\sum_{i=1}^{k}\mu_{0}^{i}-\frac{c}{k}$. When
$\nu_{i}>0$, we have $w^{*}_{i}=r$. For $\eta_{i}=0$ and $\nu_{i}=0$, we have
$0\leq w_{i}\leq r$, but this condition is not feasible as it requires the
exact identity $\mu_{0}^{i}=\lambda^{*}$. For $\eta_{i}>0$, we have
$w^{*}_{i}=0$ and $\nu_{i}=0$. Eq. 65 gives $k$ non-zeros weights which
satisfy $\mu^{i}_{0}\geq\lambda^{*}$. Therefore, we have non-zero uniform
weights $w^{i}=r$ for $i\in[1,k]$. We then apply the normalization condition
$\sum_{i=1}^{N}w_{i}=1$, which yields a uniform allocation solution with
$r_{i}^{*}=\frac{1}{k}$ for $i\in[1,k]$. It is evident that $k=N$ for a
sufficiently large $c$ such that $\sum_{i=1}^{N}(\mu_{0}^{i}-\mu_{0}^{N})<c$,
and it results in the uniform allocation $\frac{1}{N}$.

### Appendix B: Worst-case optimization of MM portfolio

The objective for the worst-case optimization of MM portfolio is given by the
following equation:

$\boldsymbol{w}^{*}=\arg\min_{\boldsymbol{w}}\left[r+\frac{b}{2}\cdot\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}\right],\,\,s.t.\,w_{i}\geq
0,\sum_{i=1}^{N}w_{i}=1,w_{i}\leq r$ (66)

The Lagrangian function is given by:

$L(w,\lambda,\eta,\nu)=r+\frac{b}{2}\cdot\boldsymbol{w}^{T}\boldsymbol{\Sigma}\boldsymbol{w}+\lambda\cdot(1-\sum_{i=1}^{N}w_{i})-\sum_{i=1}^{N}\eta_{i}\cdot
w_{i}+\sum_{i=1}^{N}\nu_{i}\cdot(w_{i}-r)$ (67)

KKT conditions:

  1. 1.

Stationarity:
$\nabla_{\boldsymbol{w}}L=b\boldsymbol{\Sigma}\boldsymbol{w}-\lambda\cdot\mathbf{1}-\boldsymbol{\eta}+\boldsymbol{\nu}=0$,
$\partial_{r}L=1-\mathbf{1}\nu=0$

  2. 2.

Primal feasibility: $w_{i}\geq 0,\sum_{i=1}^{N}w_{i}=1,w_{i}\leq c$

  3. 3.

Dual feasibility: $\eta_{i}\geq 0,\nu_{i}\geq 0$

  4. 4.

Complementary slackness: $\eta_{i}\cdot w_{i}=0,\nu_{i}\cdot(w_{i}-r)=0$

From stationarity condition, we have:

$\boldsymbol{w}=\frac{\boldsymbol{\Sigma}^{-1}}{b}\left(\lambda\cdot\mathbf{1}+\boldsymbol{\eta}-\boldsymbol{\nu}\right)$
(68)

Complementary slackness condition has three scenarios:

$1.\,\,w_{i}=0:\quad\eta_{i}\geq 0\text{ and }\nu_{i}=0$

In this scenario, the stationarity condition implies:
$(\boldsymbol{\Sigma}^{-1}(\lambda\mathbf{1}+\boldsymbol{\eta}))_{i}=0$. Let’s
denote the i-th row of the inverse of the matrix $\boldsymbol{\Sigma}$ as
$\boldsymbol{\Omega}_{i}=(\boldsymbol{\Sigma}^{-1})_{i}$. Then the condition
is $\boldsymbol{\Omega}_{i}(\lambda\mathbf{1}+\boldsymbol{\eta})=0$.

$2.\,\,0<w_{i}<r:\quad\eta_{i}=0\text{ and }\nu_{i}=0$

The stationarity condition simplifies to
$w_{i}=\frac{\lambda}{b}\left(\boldsymbol{\Omega}_{i}\cdot\mathbf{1}\right)$.

$3.\,\,w_{i}=r:\quad\eta_{i}=0\text{ and }\nu_{1i}\geq 0$

The stationarity condition implies that
$\boldsymbol{\Omega}_{i}(\lambda\mathbf{1}-\boldsymbol{\nu})=rb$

It is instructive to analyze the optimal allocation for a diagonal matrix
$\boldsymbol{\Sigma}=diag(\sigma_{1}^{2},\sigma_{2}^{2},\ldots,\sigma_{N}^{2})$.
Assume the assets are ordered in the order of increasing volatility
$\sigma_{i}$: $\sigma_{1}\leq\sigma_{2},\ldots,\leq\sigma_{N}$. According to
the slackness condition, we have three groups of indexes: $i0$ corresponds to
$w_{i0}=0$, $i+$ corresponds to $0<w_{i+}<r$ and $ir$ corresponds to
$w_{ir}=r$. Thus, we have the following condition for $\boldsymbol{\nu}$ for
each group of indexes:

$\nu_{ir}=\lambda-rb\sigma_{ir}^{2}\text{ and }\nu_{ir}\geq
0;\,\,\nu_{i0}=0;\,\,\nu_{i+}=0$ (69)

The number $k$ of weights $w_{i}=r$ is given by the maximum $k$ which still
satisfies $\lambda^{*}\geq rb\sigma_{k}^{2}$. The second stationarity
condition $\sum_{i=1}^{N}\nu_{i}=1$ is $\sum_{i=1}^{N}(\lambda-
rb\sigma_{i}^{2})^{+}=1$, where $(x)^{+}=\max(0,x)$. Thus,
$\lambda^{*}=\frac{1}{k}(1+br^{*}\sum_{i=1}^{k}\sigma_{i}^{2})$. The
conditions for $\boldsymbol{\eta}$ are:

$\eta_{ir}=0;\,,\,\eta_{i0}=-\lambda\text{ and }\eta_{i0}\geq
0;\,,\,\eta_{i+}=0$ (70)

For diagonal matrix $\boldsymbol{\Sigma}$, $\lambda^{*}$ is positive, i0 set
is empty, and $\eta_{i}=0$ for all $i$. The normalization condition
$\sum_{i=1}^{N}w_{i}=1$ implies
$kr^{*}+\lambda^{*}\sum_{i=k+1}^{N}\frac{1}{b}\frac{1}{\sigma_{i}^{2}}=1$.
Solving the above equations, we obtain the optimal values of $\lambda^{*}$,
$r^{*}$ and $k$. The optimal weights are given by Eq. 68. There are $N-k$ $i+$
weights which are inversely proportional to the variance of assets normalized
by $\lambda^{*}$ $w^{*}_{i+}=\frac{\lambda^{*}}{b}\frac{1}{\sigma_{i+}^{2}}$,
and $k$ $ir$ weights $w^{*}_{ir}=r^{*}$. The KKT-based solution impeccably
aligns with the numerical resolution computed using CVXPY.

For general $\boldsymbol{\Sigma}$, we have three group of weights: $w_{i0}=0$,
$w_{ir}=r$ and $0<w_{i+}<r$
$\boldsymbol{w}_{i+}=\frac{\lambda*}{b}\left(\boldsymbol{\Omega}_{i+}\cdot\mathbf{1}\right)$
where the weights $\boldsymbol{w}$ are a function of $\boldsymbol{\Sigma}$ and
$b$ and are obtained via the solution of the KTT equations or direct solution
obtained with the CVXPY solver.

## Appendix C: Block constraints on a sample portfolio

We examine a portfolio comprising Exxon Mobil (XOM) and Chevron (CVX) energy
stocks, as well as Pfizer (PFE) and Merck (MRK) healthcare stocks, utilizing
daily returns for the year 2022. Assuming conditional independence between the
energy and healthcare sectors, the structural information is encoded by an
adjacency matrix $\boldsymbol{A}$, containing zeros and ones:

$\boldsymbol{A}=\left(\begin{tabular}[]{rrrr}1&1&0&0\\\ 1&1&0&0\\\ 0&0&1&1\\\
0&0&1&1\\\ \end{tabular}\right)$

We estimate the precision matrix $\boldsymbol{\Theta}_{\text{sample}}$ by
inverting the sample covariance matrix
$\boldsymbol{\Theta}_{\text{sample}}=\boldsymbol{\Sigma}_{\text{sample}}^{-1}$
and the model estimation of the precision matrix
$\boldsymbol{\Theta}_{\text{model}}$ as defined in section 7 with adjacency
matrix $\boldsymbol{A}$ as follows:

$\boldsymbol{\Theta}_{sample}=\left(\begin{tabular}[]{rrrr}8974&-8400&198&-311\\\
-8400&10335&-684&-160\\\ 198&-684&5075&-3614\\\ -311&-160&-3614&9188\\\
\end{tabular}\right),\,\boldsymbol{\Theta}_{model}=\left(\begin{tabular}[]{rrrr}8438&-7827&0.0&0.0\\\
-7827&9585&0.0&0.0\\\ 0&0&4898&-3581\\\ 0&0&-3581&8962\\\
\end{tabular}\right)$

The sample covariance $\boldsymbol{\Sigma}_{sample}$ and model covariance
$\boldsymbol{\Sigma}_{model}=\boldsymbol{\Theta}_{model}^{-1}$ matrices
(multiplied by 10,000 for readability) are as follows:

$\boldsymbol{\Sigma}_{sample}=\left(\begin{tabular}[]{rrrr}4.9&4.0&0.7&0.5\\\
4.0&4.3&0.8&0.5\\\ 0.7&0.8&2.9&1.2\\\ 0.5&0.5&1.2&1.6\\\
\end{tabular}\right),\,\,\boldsymbol{\Sigma}_{model}=\left(\begin{tabular}[]{rrrr}4.9&4.0&0.0&0.0\\\
4.0&4.3&0.0&0.0\\\ 0.0&0.0&2.9&1.2\\\ 0.0&0.0&1.2&1.6\\\ \end{tabular}\right)$

The sample correlation $Corr_{sample}$ and model correlation $Corr_{model}$
matrices are as follows:

$Corr_{sample}=\left(\begin{tabular}[]{rrrr}1.00&0.88&0.19&0.19\\\
0.88&1.00&0.23&0.20\\\ 0.19&0.23&1.00&0.55\\\ 0.19&0.20&0.55&1.00\\\
\end{tabular}\right),\,\,Corr_{model}=\left(\begin{tabular}[]{rrrr}1.00&0.87&0.00&0.00\\\
0.87&1.00&0.00&0.00\\\ 0.00&0.00&1.00&0.54\\\ 0.00&0.00&0.54&1.00\\\
\end{tabular}\right)$

The sample partial correlation $\boldsymbol{\varrho}_{sample}$ and the model
partial correlation $\boldsymbol{\varrho}_{model}$ matrices are as follows:

$\boldsymbol{\varrho}_{sample}=\left(\begin{tabular}[]{rrrr}1.00&0.87&-0.03&0.03\\\
0.87&1.00&0.09&0.02\\\ -0.03&0.09&1.00&0.53\\\ 0.03&0.02&0.53&1.00\\\
\end{tabular}\right),\,\boldsymbol{\varrho}_{model}=\left(\begin{tabular}[]{rrrr}1.00&0.88&0.00&0.00\\\
0.88&1.00&0.00&0.00\\\ 0.00&0.00&1.00&0.55\\\ 0.00&0.00&0.55&1.00\\\
\end{tabular}\right)$

The off-diagonal elements of $\boldsymbol{\varrho}_{\text{sample}}$ are
smaller compared to the off-diagonal elements of $Corr_{\text{sample}}$. This
finding is consistent with our hypothesis that incorporating sector
information results in conditionally independent blocks within the portfolio.

### Appendix D: Regression for the ALD parameters

(a) The empirical relation between the asymmetric Laplace distribution
parameters for daily (blue), weekly (orange), and monthly (green) returns is
illustrated. The term $\frac{\mu}{\sigma}$ represents the location-to-scale
ratio, while $\kappa-1$ represents the skewness of the distribution. Data
points are fitted using a linear regression model, and the shaded region
indicates the error bars. The quality of fit is assessed using the $R^{2}$
score, which is presented in the legend. It is noteworthy that the slope
remains consistent across all timeframes, while the intercept demonstrates a
monotonic increase.


# RISC-V processor enhanced with a dynamic micro-decoder unit

Juliette POTTIER†, Thomas NIEDDU‡, Bertrand LE GAL‡, Sébastien PILLEMENT†, and
Maria MENDEZ REAL  †Nantes Université, CNRS, IETR-UMR 6164, F-44000 Nantes,
France
‡Université de Rennes, IRISA/INRIA lab., CNRS UMR 6074, 35042 Rennes, France

###### Abstract

For years, the open-source RISC-V instruction set has been driving innovation
in processor design, spanning from high-end cores to low-cost or low-power
cores. After a decade of evolution, RISC architectures are now as mature as
the CISC architectures popularized by industry giant Intel. Security and
energy efficiency are now joining execution speed among the design
constraints. In this article, we assess the benefits and costs associated with
integrating a micro-decoding unit inspired by CISC processors into a RISC-V
core. This unit, added in a specific pipeline stage, should enable dynamic
custom instruction sequences execution whose usage could be, for instance to
compress binaries, obfuscate behavior, etc.

###### Index Terms:

Architecture, RISC-V, ISA, CISC, micro-decoding

##  I Introduction

The field of Internet of Things (IoT) [1] and Edge Computing [2] requires the
production of a wide range of processor architectures that meet various
application constraints (computing power, cost, energy, security, etc.).
Currently, Complex Instruction Set Architecture (CISC) and Reduced Instruction
Set Architecture (RISC) processors provide high levels of performance, but in
two different ways [3]. The efficiency of CISC architectures comes from their
ability to implement complex operations with a single assembler instruction.
They provide the compiler with many macro-instructions that are efficiently
decoded and executed by the hardware architecture. These macro-instructions
offer advantages from an application perspective; for example, they allow for
optimal utilization of the hardware architecture on which the application is
executed, and they also enable a high density of binary programs. These
architectures are mainly developed by large industrial groups such as Intel or
AMD and are therefore closed. Efforts from the scientific community have
enabled the extraction of microcode from AMD K8 and K10 processors and
understanding the operation of their micro-decoding unit [4]. While being
transparent to the end user, the microcode (stored in dedicated a memory) can
be updated by the manufacturer to fix potential bugs and/or counter certain
post-fabrication attacks. This is notably what allowed Intel to address the
security vulnerabilities SPECTRE [5] and MELTDOWN [6].

RISC architectures are designed to integrate few instructions, thus lacking
the micro-decoding mechanism. The open-source RISC-V ISA (Instruction Set
Architecture), in its RVI configuration, provides the compiler with about
fifty elementary instructions [7]. Many architectures, both low-cost [8, 9,
10, 11, 12] and high-end [13, 14, 15], implement this instruction set.
Currently, high-end RISC-V processors feature 64-bit data paths, deep
pipelines, and are capable of running a Linux-type operating system thanks to
their advanced architectural optimizations. When using all these
architectures, it is up to the compiler to identify the appropriate
instruction combinations to generate efficient code. This inevitably leads to
the production of larger programs compared to their CISC counterparts, despite
the use of compressed instruction formats [16, 17]. Furthermore, unlike CISC
architectures, hardware implementations of this reduced instruction set do not
allow for post-fabrication architecture updates necessary, for example, to
address security vulnerabilities, fix bugs, better adapt to the application
domains or support new instruction set extensions.

In this paper, we investigate the introduction of a micro-decoder unit into
RISC-V processor architectures to provide flexibility. This mechanism,
introduced into the architecture, enables the decoding of custom instructions
to RISC-V instructions, or from RISC-V instructions to RISC-V instructions.
This work would allow the designer to: (1) dynamically alter the behavior of
the micro-architecture without modifying the original binary, (2) compress the
size of programs, or (3) modify the execution of an instruction. To our
knowledge, there is no existing work in the literature addressing the
integration of such a unit into a RISC processor.

The rest of the article is structured as follows. Section II presents the
study’s context and the selected RISC-V core architecture. Section III
discusses the strategies for integrating the micro-decode unit, our chosen
approach, and the necessary architectural modifications. The testing
environment, experimental methodology, and initial results are detailed in
Section IV. Section V discusses these results before concluding.

##  II Micro-decoding unit proposal

Figure 1: CV64A6 core architecture enhanced with our micro-decoding unit

Various high-end processor architectures compatible with the RISC-V
instruction set are currently available. The most popular processor cores are
BOOM [14], CVA6 [13], SHAKTI and ROCKET cores. To achieve high-performance
levels, they incorporate many pipeline stages that partition the computations
to be executed. The number of pipeline stages generally varies from 5 to a
dozen, depending on the core setup. In the most traditional configuration,
architectures supporting compressed instructions, such as the CVA6 [13],
feature 6 pipeline stages. The first 2 stages, called Frontend, handle PC
generation, instruction cache access, branch prediction and an instruction
queue. The following stage (ID) realigns compressed instructions and
translates instructions into control signals. The next stage, Issue,
integrates the scoreboard. The scoreboard controls the next 2 stages,
scheduling the operations to be executed according to available resources.
This is where the register queues (GPR and FPR) are accessed. The execution
stage (EX) takes care of the execution of the operations by the functional
units. Finally, the last stage, Commit, removes the instructions from the
pipeline and retrieves the pending results.

The addition of a CISC-like micro-decoding unit can be introduced at different
levels. The simplest strategy is to integrate it ahead of the instruction
decoder (after the fetch stage). The main advantage of this approach is that
macro-instructions can be easily detected, as they are still compressed to 32
bits. In addition, it is possible to reuse the traditional instruction decoder
to generate the control signals that will drive the architecture. However,
this approach has some notable drawbacks. The main limitation is the
impossibility of addressing additional registers outside GPR space to store
temporary computation data, thus limiting the interest of the approach.
Indeed, it discards the possibility to store temporary values. Consequently,
it was decided to add a micro-decoding unit between the Decode stage and the
Issue one, as shown in Figure 1. A specific pipeline stage, dedicated to this
task, is added to provide both (1) flexibility and modularity and (2) to avoid
impacting the critical path of the core.

##  III Microdecoding unit design

The micro-decoding unit must be able to generate sequences of $\mathbf{N_{P}}$
instructions for each of the $\mathbf{P}$ macro-instructions it needs to
handle. Different strategies are possible for generating such sequences, such
as the use of hardwired logic. However, hardware complexity will quickly
increase with the values of $\mathbf{N_{P}}$ and $\mathbf{P}$, and the
architecture of the micro-decoding stage will need to be revised with each
modification of the sequences. Moreover, this approach disables post-design
microcode upgrade. The most generic approach is based on the use of memory,
which is responsible for storing the $\mathbf{P}$ sequences of
$\mathbf{N_{P}}$ micro-instructions.

Assuming that each macro-instruction is made up of a maximum of
$\mathbf{N_{P}}$ micro-instructions, and that for address simplicity, macro-
instructions reserve $\mathbf{P}$ slots for all macro-instructions, then
memory is made up of $\mathbf{N_{P}}\times\mathbf{P}$ elements. Instructions
relating to the macro-instruction $\mathbf{idx}$ are stored in the address
range $[\mathbf{idx}\times\mathbf{N},(\mathbf{idx}+1)\times\mathbf{N}-1]$.
Each binary word in this memory represents an instruction in the sequence
$\mathbf{N_{P}}$. To reduce the hardware complexity of the microdecoding unit,
it is important that these binary words are as compact as possible. However,
they must first indicate the functional unit on which the operation is to be
performed (e.g. _ALU_ (Arithmetic-Logic Unit), _LSU_ (_Load-Store Unit_), …)
and the operation to execute (e.g. ADD, XORI, …). In the case of the CVA core,
this information is encoded using the _fu_ and _opcode_ fields, as shown in
Figure 2. Initially, we assume that the sequences of arithmetic and logical
operations for the two fields can be stored on 4 bits and 8 bits,
respectively.

Figure 2: Micro-decoding unit ROM content for macro-instruction encoding for
S-Box computation

In addition, it is necessary to memorize the source (rs1 and rs2) and
destination (rd) registers to be used by the micro-instructions. When
microdecoding a macro-instruction, the ${rd}$ register provisioned by the
compiler for macro-instruction execution may not be enough to store temporary
data. Consequently, the addresses of the rd, rs1 and rs2 registers are encoded
on 3 bits, enabling each instruction to access rs1, rs2, rd as well as 5
additional registers (added to the GPR) which are only accessible via the
microdecoder. Finally, to make the microdecoder more flexible, two additional
fields have been added. The first (imm), allows manipulation of immediate
values, while the second (skip) indicates the current instruction is the last
valid instruction in the sequence.

Without further optimization of the microdecoder’s memory structure, and
assuming that immediate data is encoded on $10$ bits, the $\mathbf{M}$ size of
each microinstruction is $32$ bits. The overall cost of memory would evolve
linearly with $\mathbf{P}$:
$\mathbf{P}\times\mathbf{N}\times\mathbf{M}=\mathbf{P}\times\mathbf{N}\times
32$ bits. This format is consistent with the size of the memory blocks
available on FPGA circuits. It should be noted that it would be possible to
implement special coding of micro-instructions (type-R, type-I, etc.) to
minimize memory complexity, but this would require decoding logic.

The microinstruction memory is controlled and interfaced with the processor
pipeline by a state machine (FSM: Finite State Machine). This machine has two
modes of operation: (1) a bypass mode that transmits conventional RISC-V
instructions to the next pipeline stage; and (2) a mode for injecting micro-
instructions into the pipeline, which is triggered as soon as a macro-
instruction is detected. The microdecoding mode includes address calculation
based on the sequence to be injected. Synchronization with other stages can be
achieved using a register barrier or even a register queue of the FIFO (First
In First Out) type.

As depicted in figure 1, the additional pipeline stage integrates 3 elements:

  1. 1.

a FIFO in charge of synchronization (acknowledgement signals) with the Issue
stage. It enables instructions to be delayed, as the microdecoder can generate
one instruction per cycle, thus saturating the scoreboard. Its depth has been
set at 2, while its width is linked to the control signals carried in the
pipeline, i.e. 364 bits. This is a FIFO of the FWFT type (First-Word-Fall-
Through), minimizing the data access penalty.

  2. 2.

an FSM associated with an address counter. In its initial state, it behaves
like a register barrier: it stores the relevant information (PC, immediate
value, rd, rs1, rs2, …) of decoded instructions. When a decoded macro-
instruction is presented, it triggers a microdecoding sequence. Based on its
$\mathbf{idx}\in\mathbf{P}$, the FSM calculates the address corresponding to
the correct sequence in memory and enters the microdecoding state. When the
sequence run is complete, the FSM returns to the (bypass) state.

  3. 3.

the memory containing the sequences of micro-instructions described in the
previous section. An example of memory from one of the test applications (AES
encryption/decryption) is shown in Figure 2. This calculation generates the
contents of the S-Box. It requires logical operations: SHIFT, AND, OR, XOR in
order to rotate on 8 bits ($\rcirclearrowleft$). This operation is broken down
into 18 RISC-V instructions. In the sequence, the registers marked $t$ refer
to temporary registers accessible only via the microdecoding unit.

##  IV Experimental results

The impact of inserting the micro-decoding unit inside a new pipeline stage
was evaluated in terms of execution time penalties of a set of test
applications. To achieve this, we modified a CV64A6 [13] core. The
architecture was simulated at cycle level using the Verilator tool (version
4.110) to accurately evaluate the execution times of various benchmarks from
the literature [18, 19, 20, 21]. Benchmark execution times were calculated
from CSR registers. The toolchain used to generate the test programs is based
on GCC compiler version 12.1.0.

The results of the cycle accurate simulations show that when the micro-
decoding unit is included in the CVA6 core, it has a slight impact on program
execution time, even if the program does not use the micro-decoding unit. A
slight variation of up to +0.3% was observed on 25 benchmarked applications,
This penalty is due to the time taken to re-fill the pipeline in the event of
incorrect branch prediction at runtime. The same protocol was followed for the
two pedagogical examples, including the S-Box generation and an AES encryption
using the microdecoding unit. In simulation, without using the microdecoder,
the blank CV64A6 completes the execution in 55368 clock cycles. On the
proposed architecture, when the micro-decoding unit is used to deliver the
18-instruction sequence, the timing is 54796 clock cycles. Performance is
better (-1% less clock cycles) when it is used. This is due to the fact that
ROM injects one instruction per clock cycle into the pipeline, whereas the
blank architecture has to fetch them from its instruction cache or even its
main memory, incurring penalties. This finding is similar when the study is
carried out directly on an FPGA board. However, unlike simulation results
cycle-accurate, those obtained on circuit FPGA are more difficult to obtain
and interpret. This is due to the OS, and more particularly the scheduler,
which due to its regular execution noises the measurements, making it
impossible to obtain accurate data when DDR memory and caches are no longer
simulated.

The impact on the core’s hardware complexity was also investigated. The
architecture was synthesized and implemented on a Kintex-7 FPGA
(XC7K325T-2FFG900C) using Xilinx’s Vivado 2020.2 tool. A summary of the
results obtained after placement and routing is provided in Table I. In this
assessment, the maximum operating frequency set at 50 MHz is not impacted by
the micro-decoding stage, demonstrating the absence of a critical path in our
stage.

| LUT | SRL | FF | BRAM36 | BRAM18 | DSP
---|---|---|---|---|---|---
baseline | 44148 | 0 | 24076 | 36 | 0 | 27
With $\mu$dec. | 45880 | 0 | 25078 | 36 | 0 | 27
(2 instr.) | +3,9% | - | +4,2% | - | - | -
With $\mu$dec. | 45922 | 0 | 25046 | 37 | 0 | 27
(32 instr.) | +3,8% | - | +3,8% | +2,6% | - | -
With $\mu$dec. | 45937 | 0 | 25068 | 38 | 0 | 27
(64 instr.) | +4,1% | - | +4,1% | +5,6% | - | -
TABLE I: Hardware complexity comparison

The first experiment illustrates the case where the micro-decoding unit has
only 2 micro-instructions. In this context, there was a slight increase in
hardware complexity of 3.9% for LUTs and 4.2% for FFs. These increases are
mainly due to the addition of the FSM, 5 registers in the GPR queue (and its
multiplexing logic) and the insertion of a FIFO of depth 2 and width 364 bits
between the microdecoding stage and the execution stage. Note that in its
version with more than 2 micro-instructions, the memory of our stage is
implemented by distributed memory. The type of implementation changes as the
size of the memory containing the micro-instructions increases. For example,
when 64 macro-instructions are stored, 2 RAM blocks are added to the LUTs and
FFs.

##  V Conclusion

In this article, a CISC-type microdecoding unit has been implemented in a
RISC-V ecosystem processor. This solution, involving the addition of an extra
pipeline slice, results in only a small latency overhead for applications not
using this functionality. Moreover, the additional costs in terms of hardware
complexity are low in comparison with the complexity of the CVA6 core. This
work shows that it is possible to hybridize a RISC-V architecture by adding a
microdecoder. This enhancement opens the way to new research perspectives for:
(1) binary size reduction, identifying recurrent instruction patterns and
incorporating them into the micro-decoding memory (2) for static code
obfuscation purposes, adding specialized instructions, akin to “black boxes,”
(3) to address security concerns for instance by injecting phantom
instructions [22, 23] into the processor pipeline and/or alter the way
calculations are performed to complicate attacks via side-channels.

## Acknowledgment

We thank the ANR which supports this work under the convention
ANR-21-CE-39-0017.

## References

  * [1] M. Hasan, “State of IOT 2022: Number of connected IOT devices growing 18% to 14.4 billion globally,” https://iot-analytics.com/number-connected-iot-devices, 2022, consulté: 29 Mars 2023.
  * [2] V. N. Chander and K. Varghese, “A soft RISC-V vector processor for Edge-AI,” in _Proceedings of VLSID_ , 2022.
  * [3] M. Burrell, _Fundamentals of Computer Architecture_. Red Globe Press London, 2004.
  * [4] P. Koppe, B. Kollenda, and al., “Reverse Engineering x86 Processor Microcode,” _CoRR_ , vol. abs/1910.00948, 2019.
  * [5] P. Kocher, J. Horn, and al., “Spectre Attacks: Exploiting Speculative Execution,” in _Proceedings of S &P’19_, 2019.
  * [6] M. Lipp, M. Schwarz, and al., “Meltdown: Reading Kernel Memory from User Space,” in _Proceedings of USENIX Security_ , 2018.
  * [7] S. L. Harris and D. Harris, _Digital Design and Computer Architecture, RISC-V Edition_. Elsevier Inc., 2021.
  * [8] Y.-H. Cheng, L.-B. Huang, and al., “RV16: An Ultra-Low-Cost Embedded RISC-V Processor Core,” _Journal of Computer Science and Technology_ , vol. 37, no. 6, pp. 1307–1319, 2022.
  * [9] K. Patsidis, D. Konstantinou, and al., “A low-cost synthesizable RISC-V dual-issue processor core leveraging the compressed Instruction Set Extension,” _Microprocess. Microsystems_ , vol. 61, pp. 1–10, 2018.
  * [10] D. A. Santos, L. M. Luza, and al., “A Low-Cost Fault-Tolerant RISC-V Processor for Space Systems,” in _Proceedings of DTIS_ , 2020.
  * [11] R. Serrano and al., “A Low-Power Low-Area SoC based in RISC-V Processor for IoT Applications,” in _Proceedings of ISOCC_ , 2021.
  * [12] N. Albartus, C. Nasenberg, and al., “On the design and misuse of microcoded (embedded) processors — a cautionary note,” in _Proceedings of USENIX Security_ , 2021.
  * [13] F. Zaruba and L. Benini, “The Cost of Application-Class Processing: Energy and Performance Analysis of a Linux-Ready 1.7-GHz 64-Bit RISC-V Core in 22-nm FDSOI Technology,” _IEEE Transactions on Very Large Scale Integration (VLSI) Systems_ , 2019.
  * [14] J. Zhao, B. Korpan, and al., “SonicBOOM: The 3rd Generation Berkeley Out-of-Order Machine,” in _Proceedings of the Workshop on Computer Architecture Research with RISC-V_ , 2020.
  * [15] D. Petrisko, F. Gilani, M. Wyse, and al., “BlackParrot: An Agile Open-Source RISC-V Multicore for Accelerator SoCs,” _IEEE Micro_ , vol. 40, no. 4, pp. 93–102, 2020.
  * [16] H. Lekatsas and W. Wolf, “Code compression for embedded systems,” in _Proceedings of DAC_ , 1998.
  * [17] H. Lozano and M. Ito, “Increasing the code density of embedded risc applications,” in _Proceedings of ISORC_ , 2016.
  * [18] B. Le Gal and C. Jego, “Softcore Processor Optimization According to Real-Applicaion Requirements,” _IEEE Embedded Systems Letters_ , vol. 5, no. 1, pp. 4–7, 2013. [Online]. Available: https://hal.science/hal-00945635
  * [19] M. Guthaus, J. Ringenberg, and al., “MiBench: A free, commercially representative embedded benchmark suite,” in _Proceedings of WWC-4 (Cat. No.01EX538)_ , 2001.
  * [20] Fossati, Luca, “TRAP benchmark suite for SystemC and TLM based Instruction Set Simulators (ISS),” http://code.google.com/p/trap-gen.
  * [21] S. Lin and D. J. Costello, _Error control coding: fundamentals and applications_. Upper Saddle River, NJ: Pearson/Prentice Hall, 2004.
  * [22] G. Leplus, O. Savry, and L. Bossuet, “Insertion of random delay with context-aware dummy instructions generator in a RISC-V processor,” in _Proceedings of HOST_ , 2022.
  * [23] M. Taram, A. Venkat, and D. Tullsen, “Mitigating Speculative Execution Attacks via Context-Sensitive Fencing,” _IEEE Design & Test_, vol. 39, no. 4, pp. 49–57, 2022.