1. Introduction

Portfolio selection aims at either maximizing the return or minimizing the risk. In 1952, Markowitz (1952) suggests to select the portfolio by minimizing the standard deviation at a given expected return under the assumption that asset returns are normally distributed. This means that standard deviation is chosen as the risk measure. Markowitz's work laid down the cornerstone for modern portfolio selection theory framework.

Risk measures and probability distributions are two important constituents of the portfolio selection theory. Traditional Markowitz's model (Markowitz, 1952) is established based on normality assumption and standard deviation is chosen as the risk measure.

One disadvantage of taking standard deviation (StD) as a risk measure is that the loss in the extreme cases tends to be underestimated. To overcome such a difficulty, the idea of Value at Risk (VaR) is also widely used in practice. Artzner et al. (1999) suggests that desirable risk measure should be "coherent". However, VaR does not fulfill the subadditivity condition as required by the definition of "coherence". Yiu (2004) proposed an optimal portfolio selection under Value-at-Risk. On the other hand, Expected Shortfall (ES) is coherent as a popular risk measure for portfolio selection that aims at averaging the tail uncertainties.

It is well-known that financial data cannot be described satisfactorily by normal distribution. The normality assumption is restrictive and is generally violated due to financial market uncertainties and managers' risk aversion. As Behr and Ptter (2009) pointed out, alternatives for multivariate normal distribution are necessary for portfolio selection. A desirable alternative model should be able to explain tail heaviness, skewness, and excess kurtosis. Various heavy tailed distributions have been applied to portfolio selection problems. Among these, Mandelbrot (1997) concluded that the daily rate of return of stock price data exhibit heavy tailed distributions; Hu and Kercheval (2010) apply multivariate skewed $t$ and student $t$ distribution for efficient frontier analysis; Generalized hyperbolic distribution is extensively studied in (Behr and Ptter, 2009; Eberlein, 2001; Hellmich and Kassberger, 2011; Hu and Kercheval, 2007; Surya and Kurniawan, 2014; Socgnia and Wilcox, 2014), with special cases including hyperbolic distribution (Bingham and Kiesel, 2001; Eberlein and Keller, 1995), Variance Gamma distribution (Seneta, 2004), Normal Inverse Gaussian distribution (Barndor-Nielsen, 1995), etc.

Recently, Asymmetric Laplace distribution has received various attention in the literature, to name a few, (Ayebo and Kozubowski, 2003; Kollo and Srivastava, 2005; Kozubowski and Podgrski, 1999; Kozubowski and Podgrski, 2001; Punathumparambath, 2012). Compared to Normal distribution, the Asymmetric Laplace distribution describes asymmetry, steep peak, and tail heaviness better. Portfolio selection models are extensively studied under Asymmetric Laplace framework. Zhu (2007), Kozubowski and Podgrski (2001) apply Asymmetric Laplace distribution to financial data. By assuming that the asset data is generated from autoregressive moving average (ARMA) time series models with Asymmetric Laplace noise, Zhu (2007) establish the asymptotic inference theory under very mild conditions and present methods of computing conditional Value at Risk (CVaR). Zhao et al. (2015) further propose a so-called mean-CVaR-skewness portfolio selection strategy under Asymmetric Laplace distribution, this model can be further transformed to quadratic programming problem with explicit solutions.

In this paper, we extended Hu's work (Hu, 2010) to Asymmetric Laplace framework. We first derived the equivalence of mean-VaR/ES/Std-skewness-kurtosis models, and show that these models can be reduced to quadratic programming problem. Since Zhao et al. (2015) utilized moment estimation for parameter estimation of Asymmetric Laplace distribution which less efficient compare to maximum likelihood estimation. Taken into consideration of the normal mean-variance mixture of Asymmetric Laplace distribution, followed by Expectation-Maximization algorithm for multivariate Laplace distribution in Arslan (2010), we derived the EM algorithm for Asymmetric Laplace Distributions that outperforms moment estimation in Zhao et al. (2015). The advantage of the proposed EM algorithm is to alleviate the complicated calculation of Bessel function. This improves many existing methods of estimating Asymmetric Laplace distributions, for example, Hrlimann (2013), Kollo and Srivastava (2005) and Visk (2009). Extensive simulation studies and efficient frontier analysis are complemented to confirm that our algorithm performs better than moment estimation for parameter estimation.

The rest of the article is organized as follows. In Section 2, properties of Asymmetric Laplace distributions and coherent risk measures are summarized. In Section 3, portfolio selection under Asymmetric Laplace framework are derived, complement with Expectation-Maximization (EM) algorithm for parameter estimation of Asymmetric Laplace distributions. In Section 4, simulation studies are provided to show the efficiency of the Expectation-Maximization procedure. Section 5 presents real data analysis of Asymmetric Laplace Distributions based portfolio selection models, followed by conclusive remarks in Section 6.

---

2. Preliminary knowledge

---

2.1. Asymmetric Laplace distribution

Kotz et al. (2001) proposed the Asymmetric Laplace Distribution with density function

$f({\boldsymbol{x}}) = \frac{ 2e^{{\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}} } { (2\pi)^{n/2}|{\bf{\Sigma}}|^{1/2} } \big{(} \frac{{\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{x}} }{ 2+{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}} \big{)} ^{v/2} K_{v} \big{(} \sqrt{ (2+{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}})({\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{x}}) } \big{)},$

(2.1)

denoted as ${\boldsymbol{X}}\sim AL_n({\boldsymbol{\mu}}, {\bf{\Sigma}})$. Here, $n$ is the dimension of random vector ${\boldsymbol{X}}$, $v = (2-n)/2$ and $K_{v}(u)$ is the modified Bessel function of the third kind with the following two popular representations:

$K_v(u) = \frac{1}{2} \big(\frac{u}{2}\big)^v\int_{0}^{\infty}t^{-v-1}\exp\Big{\{} -t-\frac{u^2}{4t} \Big{\}} dt, \quad u>0,$

(2.2)

$K_v(u) = \frac{(u/2)^v\Gamma(1/2)}{\Gamma(v+1/2)}\int_{1}^{\infty} e^{-ut}(t^2-1)^{v-1/2}dt, \quad u>0, v\geq -1/2.$

(2.3)

When ${\boldsymbol{\mu}} = {\boldsymbol{0}}_n$, we can obtain Symmetric Laplace distribution $SL\, ({\bf{\Sigma}})$ with density

$f({\boldsymbol{x}}) = 2(2\pi)^{-n/2}|{\bf{\Sigma}}|^{-1/2}\Big({\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{x}}/2 \Big{)}^{\nu/2} K_{\nu}\Big(\sqrt{ 2{\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{x}} } \Big).$

When $n = 1$, we have ${\bf{\Sigma}} = \sigma_{11} = \sigma\, .$ In such cases, (2.1) becomes the univariate Laplace distribution $AL_1(\mu, \sigma)$ distribution with parameters $\mu$ and $\sigma$. The corresponding density function is

$f(x) = \frac{1}{\gamma} \exp\Big{\{} -\frac{|x|}{\sigma^2}\big[\gamma-\mu\cdot \text{sign}(x)\big] \Big{\}} \quad \text{ with } \quad \gamma = \sqrt{\mu^2+2\sigma^2}.$

(2.4)

The symmetric case ($\mu = 0$) leads to the univariate Laplace distribution $SL_1\, (0, \sigma)$.

Figure 1 displays plot of symmetric densities and AL densities. Symmetric densities including standard normal distribution, student $t$ distribution with 2 degrees of freedom, and univariate symmetric Laplace distribution, denoted as $N\, (0, 1), t\, (2), SL_1\, (0, 1)$. The student $t$ distribution possesses heavier tail than normal distribution, whereas $SL_1\, (0, 1)$ distribution imposes greater peakedness and heavier tail than normal case. As for plots of AL densities, when $\mu>0$, the density skews to the right. On the other hand, when $\mu < 0$, the density skews to the left.

Important results of univariate and multivariate asymptotic Laplace distributions that will be used later on are presented below.

Proposition 2.1. (See Kotz, 2001)

(1). If ${\boldsymbol{X}} = (X_1, \cdots, X_n)$ follows multivariate Asymmetric Laplace distribution, i.e., ${\boldsymbol{X}}\sim AL_n\, ({\boldsymbol{\mu}}, {\bf{\Sigma}})$, $n$ is the number of securities. The linear combination ${\boldsymbol{w}}'{\boldsymbol{X}} = w_1X_1+\cdots+w_nX_n$ follows univariate Asymmetric Laplace distribution, i.e. ${\boldsymbol{w}}'{\boldsymbol{X}}\sim AL_1\, (\mu, \sigma)$, with $\mu = {\boldsymbol{w}}'{\boldsymbol{\mu}}\, , \sigma = \sqrt{{\boldsymbol{w}}'{\bf{\Sigma}}{\boldsymbol{w}}}\, , {\boldsymbol{w}} = (w_1, \cdots, w_n)'$.

(2). Assume that univariate random variable $Y \sim AL_1\, (\mu, \sigma)$. To measure the asymmetry and peakedness of the distribution, define the skewness (Skew[Y]) and kurtosis (Kurt[Y]) as the third and fourth standardized moment of a random variable $Y\, .$ Then,

$\begin{eqnarray*} \text{Skew}[Y] & = & \frac{ \mathbb{E}(Y-\mathbb{E}Y)^3 }{\big[ \mathbb{E}(Y-\mathbb{E}Y)^2 \big]^{3/2}} = \frac{2\mu^3+3\mu\sigma^2}{(\mu^2+\sigma^2)^{3/2}}, \\ {\text{Kurt}[Y]} & = & \frac{\mathbb{E}(Y-\mathbb{E}Y)^4}{\big[\text{Var}(Y)\big]^2} = \frac{ 9\mu^4+6\sigma^4+18\mu^2\sigma^2 }{ (\mu^2+\sigma^2)^2 }. \end{eqnarray*}$

(3). Let ${\boldsymbol{X}} = (X_1, X_2, \cdots, X_n)\sim AL_n\, ({\boldsymbol{\mu}}, {\bf{\Sigma}})$. Then the first and second order moments of ${\boldsymbol{X}}$ are

$\mathbb{E}({\boldsymbol{X}}) = {\boldsymbol{\mu}} \quad and \quad \text{Cov}\, ({\boldsymbol{X}}) = {\bf{\Sigma}}+{\boldsymbol{\mu}}'{\boldsymbol{\mu}}.$

(4). The Asymmetric Laplace distribution can be represented as a mixture of normal vector and a standard exponential variable, i.e., ${\boldsymbol{X}} \sim AL_n\, ({\boldsymbol{\mu}}, {\bf{\Sigma}})$ can be represented as

${\boldsymbol{X}} = {\boldsymbol{\mu}} Z+Z^{1/2}{\boldsymbol{Y}},$

where ${\boldsymbol{Y}}\sim N_n({\boldsymbol{0}}, {\bf{\Sigma}})\, , Z\sim \text{Exp}(1)$. This indicate that we can simulate multivariate Asymmetric Laplace random vector ${\boldsymbol{X}}\sim AL_n({\boldsymbol{\mu}}, {\bf{\Sigma}})$ as follows:

1. Generate a multivariate normal variable ${\boldsymbol{Y}}\sim N_n\, ({\boldsymbol{0}}, {\bf{\Sigma}})$;

2. Generate a standard exponential variable $Z\sim \text{Exp}\, (1)$;

3. Construct Asymmetric Laplace random vector as ${\boldsymbol{X}} = {\boldsymbol{\mu}} Z+Z^{1/2}{\boldsymbol{Y}}$.

Figure 2 displays several realizations of bivariate Asymmetric Laplace distribution with different levels of asymmetry and peakedness.

---

2.2. Risk measures

Since mean and covariance matrix cannot be used to characterize non-Gaussian distribution, alternative risk measures are necessary for portfolio selection problems. Artzner et al. (1999) suggests that a desirable risk measure should be defined fulfilling certain properties and such a risk measure is said to be coherent.

A risk measure $\phi$ that maps a random variable to a real number is coherent if it satisfies the following conditions:

1). Translation invariance: $\phi(l+h) = \phi(l)+h$, for all random losses $l$ and all $h\in\mathbb{R}$;

2). Subadditivity: $\phi(l+h)\leq \phi(l)+\phi(h)$, for all random losses $l, h$;

3). Positive homogeneity: $\phi(\lambda l) = \lambda \phi(l)$ for all random losses $l$ and all $\lambda>0$;

4). Monotonicity: $\phi(l_1)\leq \phi(l_2)$ for all random losses $l_1\, , l_2$ with $l_1\leq l_2$ almost surely.

Standard deviation is not coherent in general excepting the Gaussian cases. VaR is coherent when the underlying distribution is elliptically distributed. Expected Shortfall, or the so-called conditional value at risk (CVaR) is a coherent risk measure since it always satisfies subadditivity, monotonicity, positive homogeneity, and convexity. For any fixed $\alpha\in(0, 1)$, $\alpha$-$\text{VaR}$ is the $\alpha$-quantile loss while $\alpha$-ES is the average of all $\beta$-VaR for $\beta\in(\alpha, 1)$. Both VaR and CVaR measure the potential maximal loss. VaR and ES can be written as

$\text{VaR}_{\alpha} = F^{-1}(\alpha)\quad \mbox{and} \quad \text{ES}_{\alpha} = E[L|L\leq \text{-VaR}_{\alpha}] = - \frac{1}{\alpha}\int_{-\infty}^{-\text{VaR}_\alpha}\text{VaR}_{\beta}d\beta,$

where $F(\cdot)$ is the cumulative distribution function of loss $L$ and $\text{ES}_\alpha$ is the expected loss above $\text{VaR}_\alpha$. Thus, the estimation process are

$\int_{-\infty}^{-\text{VaR}_\alpha} f_X(x)dx = \alpha\quad \mbox{and} \quad \text{ES}_{\alpha} = - \frac{1}{\alpha}\int_{-\infty}^{-\text{VaR}_\alpha} xf(x)dx.$

(2.5)

Under normality assumption, VaR$_\alpha$ and ES$_\alpha$ are

$\begin{eqnarray*} && \text{VaR}_\alpha = \mu+\sigma \Phi^{-1}(1-\alpha), \\ && \text{ES}_\alpha = \mu+\sigma\frac{\psi(\Phi^{-1}(1-\alpha))}{\alpha}. \end{eqnarray*}$

where $\psi(\cdot)$ as the normal density distribution, and $\Phi^{-1}(\cdot)$ is the quantile distribution.

As shown in Hu et al. (2010), portfolio selected by minimizing standard deviation, VaR$_\alpha$, and ES$_\alpha$ are the equivalent under elliptical distribution assumption.

It is well-documented that asset securities are not normally distributed. As an alternative to Gaussian distribution, Asymmetric Laplace distribution exhibits tail-heaviness, skewness, and peakedness.

---

3. Portfolio selection under ALD framework

Let ${\boldsymbol{X}} = (X_1, X_2, \cdots, X_n)\sim AL_n({\boldsymbol{\mu}}_{\mathcal{P}}, {\bf{\Sigma}}_{\mathcal{P}})$ be the return vectors of $n$ securities, and ${\boldsymbol{w}} = (w_1, w_2, \cdots, w_n)'$ is the allocation weight vector. Then, the portfolio is defined as

$\mathcal{P}\, ({\boldsymbol{w}}) = {\boldsymbol{w}}'{\boldsymbol{X}} = \sum\limits_{i = 1}^{n}w_iX_i .$

According to Proposition 2.1 (2), $\mathcal{P}\, ({\boldsymbol{w}})\sim AL_1\, (\mu, \sigma)$ with $\mu = {\boldsymbol{w}}'{\boldsymbol{\mu}}\, , \sigma = \sqrt{{\boldsymbol{w}}'{\bf{\Sigma}}{\boldsymbol{w}}}$.

From Theorem 3.1–3.2 below, in order to select a portfolio under Asymmetric Laplace distribution, it suffices to obtain the unknown parameters ${\boldsymbol{\mu}}_{\mathcal{P}}$ and ${\bf{\Sigma}}_{\mathcal{P}}$. Thus portfolio selection models under Asymmetric Laplace distribution lead to parameter estimation for $AL_n({\boldsymbol{\mu}}_{\mathcal{P}}, {\bf{\Sigma}}_{\mathcal{P}})$. Zhao et al. (2015) proposed the multi-objective portfolio selection model under Asymmetric Laplace framework and derived the simplified model that can be reformulated as quadratic programming problem. However, to estimate the unknown parameters, the authors adopt a moment estimation method that is less efficient compared to maximum likelihood method. Since Asymmetric Laplace distribution can be represented as a mixture of exponential distribution and multivariate normal distribution, we derived the Expectation-Maximization algorithm for parameter estimation of Asymmetric Laplace distribution. The algorithm for estimating these unknown parameters is discussed in Section 3.2.

---

3.1. Portfolio selection theorems

Theorem 3.1. Let ${\boldsymbol{X}} = (X_1, \cdots, X_n) \sim AL_n({\boldsymbol{\mu}}_{\mathcal{P}}, {\bf{\Sigma}}_{\mathcal{P}})$ be a $n$-dimensional random vector that follow multivariate Asymmetric Laplace distribution, each element ($X_i, i = 1, 2, \cdots, n$) represent a stock. Let ${\boldsymbol{w}}$ be the weight vector and $\mathcal{P}({\boldsymbol{w}}) = {\boldsymbol{w}}'{\boldsymbol{X}} = \sum_{i = 1}^{n}w_iX_i$ be the portfolio. Then, under Asymmetric Laplace framework, risk measures of StD, VaR$_\alpha$, and ES$_\alpha$ at $\alpha\in(0, 1)$ level formulated as

$\begin{eqnarray*} && \text{Standard Deviation: }\quad \text{StD}\, \big(\mathcal{P}({\boldsymbol{w}})\big) = \frac{\sigma}{\sqrt{2}}; \\ && \text{Value at Risk: } \quad \qquad \text{VaR}_\alpha\, \big(\mathcal{P}({\boldsymbol{w}})\big) = -\frac{\sigma^2}{\gamma+\mu}\ln \frac{\alpha\gamma(\gamma+\mu)}{\sigma^2}; \\ && \text{Expected Shortfall: }\quad \text{ES}_\alpha\, \big(\mathcal{P}({\boldsymbol{w}})\big) = \frac{\sigma^2}{\gamma+\mu}-\frac{\sigma^2}{\gamma+\mu}\ln \frac{\alpha\gamma(\gamma+\mu)}{\sigma^2}. \end{eqnarray*}$

Here, $\mu = {\boldsymbol{w}}'{\boldsymbol{\mu}}_{\mathcal{P}} = \text{mean}\, \big(\mathcal{P}({\boldsymbol{w}})\big)\, , \sigma = \sqrt{{\boldsymbol{w}}'{\bf{\Sigma}}_{\mathcal{P}}{\boldsymbol{w}}} = \text{std}\, \big(\mathcal{P}({\boldsymbol{w}})\big)$ and $\gamma = \sqrt{\mu^2+2\sigma^2}$.

Proof. Let ${\boldsymbol{\mu}}_{\mathcal{P}} = (\mu_1, \cdots, \mu_n)$ be the mean return vector of the securities $(X_1, \cdots, X_n)$ and $\Sigma_{\mathcal{P}} = \big(\sigma_{\mathcal{P}} \big)_{i, j = 1}^{p}$ be the scale matrix of $(X_1, \cdots, X_n)$. Denote the allocation vector by ${\boldsymbol{w}} = (w_1, \cdots, w_n)'$. Then, the portfolio $\mathcal{P}({\boldsymbol{w}}) = \sum_{i = 1}^{n}w_iX_i$ follows univariate Asymmetric Laplace distribution with

$\mathcal{P}({\boldsymbol{w}}) = \sum\limits_{i = 1}^{n}w_iX_i\sim AL_1(\mu, \sigma) \quad \text{ with } \quad \mu = \sum\limits_{i = 1}^{n}\mu_iw_i \, , \sigma = \Big(\sum\limits_{i = 1}^{n}\sum\limits_{j = 1}^{n}\sigma_{\mathcal{P}_{ij}}w_iw_j \Big)^{1/2}.$

If ${\boldsymbol{\mu}} = {\boldsymbol{0}}_n$, the univariate symmetric Asymmetric Laplace distribution becomes $AL_1(0, \sigma)$ with density

$g(x) = \frac{1}{\gamma}\exp\Big\{ -\frac{|x|}{\sigma^2}\gamma \Big\} \quad \text{ with } \quad \gamma = \sqrt{2}\sigma.$

Thus, standard deviation (StD) of portfolio $\mathcal{P}({\boldsymbol{w}}) = {\boldsymbol{w}}'{\boldsymbol{X}}$ is

$\text{StD}\, (\mathcal{P}({\boldsymbol{w}})) = \int_{-\infty}^{+\infty} \frac{1}{\gamma}|x|\exp\Big\{ -\frac{\gamma}{\sigma^2}|x| \Big\} dx = 2\int_{0}^{+\infty} \frac{x}{\gamma}\exp\Big\{ -\frac{\gamma}{\sigma^2}x \Big\} dx = \frac{2\sigma^4}{\gamma^3} = \frac{\sigma}{\sqrt{2}}.$

According to the definition of VaR$_\alpha$ and ES$_\alpha$ as defined in (2.5) and univariate Asymmetric Laplace density (2.4), we have

$\begin{eqnarray*} && -\int_{-\infty}^{-\text{VaR}_\alpha} \frac{1}{\gamma}\exp\Big{\{} -\frac{|x|}{\sigma^2} \big[ \gamma-\mu\cdot\text{sgn}(x)\big] \Big{\}} dx = \alpha, \\ && \frac{\sigma^2}{\gamma(\gamma+\mu)}\exp\Big{\{} -\frac{\gamma+\mu}{\sigma^2}\text{VaR}_\alpha \Big{\}} = \alpha. \\ \end{eqnarray*}$

Thus, VaR$_\alpha$ and ES$_\alpha$ are

$\begin{eqnarray*} && \text{VaR}_\alpha\, \big(\mathcal{P}({\boldsymbol{w}})\big) = -\frac{\sigma^2}{\sqrt{\mu^2+2\sigma^2}+\mu}\ln\frac{\alpha(\mu^2+2\sigma^2+\mu\sqrt{\mu^2+2\sigma^2})}{\sigma^2} = -\frac{\sigma^2}{\gamma+\mu}\ln\frac{\alpha\gamma(\gamma+\mu)}{\sigma^2}; \\ && \text{ES}_\alpha\, \big(\mathcal{P}({\boldsymbol{w}})\big) = -\frac{1}{\alpha}\int_{-\infty}^{-\text{VaR}_\alpha}xf_X(x)dx = -\frac{1}{\alpha}\int_{-\infty}^{-\text{VaR}_\alpha}x\frac{1}{\gamma}\exp\Big{\{} -\frac{|x|}{\sigma^2}\big[\gamma-\mu\cdot\text{sgn}(x)\big] \Big{\}} dx\\ && \qquad = \frac{\sigma^2}{\mu+\sqrt{\mu^2+2\sigma^2}} - \frac{\sigma^2}{\mu+\sqrt{\mu^2+2\sigma^2}}\ln\Big\{ 2\alpha+\frac{\alpha(\mu^2+\mu\sqrt{\mu^2+2\sigma^2})}{\sigma^2} \Big\} \\ && \qquad = \frac{\sigma^2}{\gamma+\mu}-\frac{\sigma^2}{\gamma+\mu}\ln\frac{\alpha\gamma(\gamma+\mu)}{\sigma^2}. \end{eqnarray*}$

Then we have the following theorem.

Theorem 3.2. Let ${\boldsymbol{X}}\sim AL_n({\boldsymbol{\mu}}_{\mathcal{P}}, {\bf{\Sigma}}_{\mathcal{P}})$. Then, portfolio $\mathcal{P}({\boldsymbol{w}}) = {\boldsymbol{w}}'{\boldsymbol{X}}$ with following models based on $ES_{\alpha}\, ,$ $VaR_{\alpha}\, ,$ and $StD$ (as defined in Theorem 3.1)

$\begin{eqnarray*}\label{opt1} && \min\limits_{{\boldsymbol{w}}}\text{ES}_{\alpha}\, \big(\mathcal{P}({\boldsymbol{w}})\big) \text{ or } \min\limits_{{\boldsymbol{w}}} \text{VaR}_{\alpha}\, \big(\mathcal{P}({\boldsymbol{w}})\big) \text{ or } \min\limits_{{\boldsymbol{w}}} \text{StD}\, \big(\mathcal{P}({\boldsymbol{w}})\big) \\ && \max\limits_{{\boldsymbol{w}}} \quad \text{Skew[} \mathcal{P}({\boldsymbol{w}}) \text{]} = \frac{2\mu^3+3\mu\sigma^2}{(\mu^2+\sigma^2)^{3/2}} \\ && \max\limits_{{\boldsymbol{w}}} \quad \text{Kurt[} \mathcal{P}({\boldsymbol{w}}) \text{]} = \frac{ 9\mu^4+6\sigma^4+18\mu^2\sigma^2 }{ (\mu^2+\sigma^2)^2 } \\ && \text{s.t.} \qquad {\boldsymbol{w}}'{\boldsymbol{\mu}} = r_0\, , {\boldsymbol{w}}'{\boldsymbol{1}} = 1. \end{eqnarray*}$

are equivalent. Here, $\mu = {\boldsymbol{w}}'{\boldsymbol{\mu}}_{\mathcal{P}} = \text{mean}\, \big[\mathcal{P}({\boldsymbol{w}})\big]\, , \sigma = \sqrt{ {\boldsymbol{w}}'{\bf{\Sigma}}_{\mathcal{P}}{\boldsymbol{w}} } = \text{std}\, \big[\mathcal{P}({\boldsymbol{w}})\big]\, , {\boldsymbol{w}} = (w_1, w_2, \cdots, w_n)'$.

Proof. Let $g(\mu, \sigma) = \frac{\sigma^2}{\mu+\sqrt{\mu^2+2\sigma^2}}$. Then, ES$_\alpha$[$\mathcal{P}({\boldsymbol{w}})$] and VaR$_\alpha$[$\mathcal{P}({\boldsymbol{w}})$] are

$\begin{eqnarray*} && \text{VaR}_\alpha[\mathcal{P}({\boldsymbol{w}})] = - g(\mu, \sigma)\ln \big(2\alpha+\frac{\alpha\mu}{g(\mu, \sigma)} \big) = -g(\mu, \sigma) \big[ \ln \alpha + \ln\big(2+\frac{\mu}{g(\mu, \sigma)} \big) \big], \\ && \text{ES}_\alpha[\mathcal{P}({\boldsymbol{w}})] = g(\mu, \sigma) - g(\mu, \sigma)\ln \big(2\alpha+ \frac{\alpha\mu}{g(\mu, \sigma)} \big) = (1-\ln\alpha)g(\mu, \sigma) - g(\mu, \sigma)\ln(2+\frac{\mu}{g(\mu, \sigma)}). \end{eqnarray*}$

Differentiating the above expressions with respect to $\sigma$, we have

$\begin{eqnarray*} && \frac{\partial\text{VaR}_\alpha[\mathcal{P}({\boldsymbol{w}})] }{\partial\sigma} = \frac{\partial g(\mu, \sigma)}{\partial\sigma}\left[ -\ln\alpha -\ln\big(2+\frac{\mu}{g(\mu, \sigma)} \big) + \frac{ \frac{\mu}{g(\mu, \sigma)} }{2+\frac{\mu}{g(\mu, \sigma) }} \right] > 0, \\ && \frac{\partial\text{ES}_\alpha[\mathcal{P}({\boldsymbol{w}})] }{\partial\sigma} = \frac{\partial g(\mu, \sigma)}{\partial\sigma} \left[ 1-\ln\alpha-\ln\big(2+\frac{\mu}{g(\mu, \sigma)} \big) + \frac{ \frac{\mu}{g(\mu, \sigma)} }{2+\frac{\mu}{2+g(\mu, \sigma)}} \right] > 0, \end{eqnarray*}$

where

$\frac{\partial g(\mu, \sigma)}{\partial \sigma} = \frac{\partial\Big[ \frac{\sigma^2}{\mu+\sqrt{\mu^2+2\sigma^2}} \Big]}{\partial\sigma} = 2\sigma \frac{ \mu+\frac{\mu^2+\sigma^2}{\sqrt{\mu^2+2\sigma^2}} }{ (\mu^2+\sqrt{\mu^2+2\sigma^2})^2 } > 0 .$

The derivative of skewness measure with respect to $\sigma$ is

$\frac{\partial \text{Skew[}\mathcal{P}({\boldsymbol{w}})\text{]}}{\partial\sigma} = \frac{\partial\Big[ \frac{2\mu^3+3\mu\sigma^2}{(\mu^2+\sigma^2)^{3/2}} \big]}{\partial\sigma} = \frac{-3\mu\sigma^3}{ (\mu^2+\sigma^2)^{5/2} } < 0.$

The derivative of kurtosis measure with respect to $\sigma$ is

$\frac{\partial \text{Kurt[}\mathcal{P}({\boldsymbol{w}})\text{]}}{\partial\sigma} = \frac{ \frac{ 9\mu^4+6\sigma^4+18\mu^2\sigma^2 }{ (\mu^2+\sigma^2)^2 } }{\partial \sigma} = \frac{ -12\mu^4\sigma^3-12\mu^2\sigma^5 }{ (\mu^2+\sigma^2)^4 } <0.$

The monotonicity of VaR$_\alpha$[$\mathcal{P}({\boldsymbol{w}})$], ES$_\alpha$[$\mathcal{P}({\boldsymbol{w}})$], Skew[$\mathcal{P}({\boldsymbol{w}})$], and Kurt[$\mathcal{P}({\boldsymbol{w}})$] with respect to $\sigma$ indicate that the portfolio selection problems based on these risk measures are equivalent. This means that minimizing VaR$_\alpha$[$\mathcal{P}({\boldsymbol{w}})$], ES$_\alpha$[$\mathcal{P}({\boldsymbol{w}})$], StD[$\mathcal{P}({\boldsymbol{w}})$] are equivalent to minimizing ${\boldsymbol{w}}'{\bf{\Sigma}}_{\mathcal{P}}{\boldsymbol{w}}$.

---

3.2. Parameter estimation of Asymmetric Laplace distribution

Assume ${\boldsymbol{X}} = (X_1, X_2, \cdots, X_n)\sim AL_n({\boldsymbol{\mu}}, {\bf{\Sigma}})$. Let ${\boldsymbol{x}}_1\, , {\boldsymbol{x}}_2\, , \cdots, {\boldsymbol{x}}_T \in \mathbb{R}^n$ be the $T$ observations. We aim at fitting a multivariate Asymmetric Laplace distribution $\mathcal{AL}_n({\boldsymbol{\mu}}, {\bf{\Sigma}})$ with unknown parameters ${\boldsymbol{\mu}}, {\bf{\Sigma}}$.

Hrlimann (2013), Kollo and Srivastava (2005), Visk (2009) consider moment matching methods that is less efficient than maximum likelihood estimation. Kotz et al. (2002) and Kotz et al. (2001) presented the maximum likelihood estimators for parameter estimation of Asymmetric Laplace distributions. However, maximum likelihood estimation require computation of complicated Bessel function. Thus we derived the expectation-maximization algorithm for parameter estimation of Asymmetric Laplace distribution.

---

3.2.1. Moment estimation

As Zhao et al. (2015) pointed out, according to Proposition 2.1 (3), Asymmetric Laplace distribution can be estimated via moment method (Moment-AL) with

$\hat{{\boldsymbol{\mu}}} = \bar{{\boldsymbol{x}}} \quad\mbox{and}\quad \hat{{\bf{\Sigma}}} = \text{cov}({\boldsymbol{X}})-\hat{{\boldsymbol{\mu}}}'\hat{{\boldsymbol{\mu}}},$

where $\bar{{\boldsymbol{x}}} = \frac{1}{n}\sum_{i = 1}^{n}{\boldsymbol{x}}_i\, , \text{Cov}({\boldsymbol{X}}) = \sum_{i = 1}^{n}({\boldsymbol{x}}_i-\bar{{\boldsymbol{x}}})^T({\boldsymbol{x}}_i-\bar{{\boldsymbol{x}}})$.

---

3.2.2. Maximum likelihood estimation

Consider sample points ${\boldsymbol{x}}_1, {\boldsymbol{x}}_2, \cdots, {\boldsymbol{x}}_n$ and density function of Asymmetric Laplace distribution as defined in (2.1). Taken logarithm with respect to likelihood function, the log-likelihood is

$\begin{eqnarray*} \ell({\boldsymbol{\mu}}, {\bf{\Sigma}}) & = & \ln L({\boldsymbol{\mu}}, {\bf{\Sigma}}) = \sum\limits_{t = 1}^{T} \ln f({\boldsymbol{x}}_t;{\boldsymbol{\mu}}, {\bf{\Sigma}}) \\ & = & \sum\limits_{t = 1}^{T}{\boldsymbol{x}}_t{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}} + T\ln 2 - \frac{Tn}{2} \ln(2\pi) - \frac{T}{2} \ln \big(|{\bf{\Sigma}}|\big) +\frac{\nu}{2} \sum\limits_{t = 1}^{T}\ln \big({\boldsymbol{x}}_t'{\bf{\Sigma}}^{-1}{\boldsymbol{x}}_t \big) - \\ && \frac{\nu T}{2} \ln \big(2+{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}} \big) + \sum\limits_{t = 1}^{T} \ln K_{v} \Big{\{} \sqrt{ (2+{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}})({\boldsymbol{x}}_t'{\bf{\Sigma}}^{-1}{\boldsymbol{x}}_t) } \Big{\}}. \end{eqnarray*}$

Generally, we can directly maximize the log-likelihood function $\ell({\boldsymbol{\mu}}, {\bf{\Sigma}})$ with respect to parameters ${\boldsymbol{\mu}}, {\bf{\Sigma}}$ and thus obtain the maximum likelihood estimator. Unfortunately, the density function involves modified Bessel function of the third kind with density (2.2), (2.3) that are too complex and complicated for numerical maximization. However, Gaussian-Exponential mixture representation of the Asymmetric Laplace distribution allows us to employ the expectation-maximization algorithm without involving modified Bessel functions.

---

3.2.3. Expectation-maximization algorithm

Then we derive the Expectation-Maximization algorithm for parameter estimation of multivariate Asymmetric Laplace Distribution (mALD), we follow the EM derivation for Multivariate Skew Laplace distribution in Arslan (2010).

Let ${\boldsymbol{X}} = (X_1, X_2, \cdots, X_n)$ be Asymmetric Laplace distributed random vector. Proposition 2.1 suggests that ${\boldsymbol{X}}$ can be generated from a latent random variable $Z = z$ through multivariate Gaussian distribution with $z{\boldsymbol{\mu}}, z{\bf{\Sigma}}$, i.e. $X|Z = z \sim N_{n}\, (z{\boldsymbol{\mu}}, z{\bf{\Sigma}})\, , Z\, \sim \text{Exp}\, (1)$ with density

$\begin{eqnarray*} f_{{\boldsymbol{X}}|Z}({\boldsymbol{x}}, z) & = & \frac{1}{(2\pi)^{n/2}|z{\bf{\Sigma}}|^{1/2}} \exp \Big{\{} -\frac{1}{2}({\boldsymbol{x}}-z{\boldsymbol{\mu}})'(z{\bf{\Sigma}})^{-1}({\boldsymbol{x}}-z{\boldsymbol{\mu}}) \Big{\}}, \\ f_{Z}(z) & = & e^{-z}1_{ \{z\geq 0\} }. \end{eqnarray*}$

Thus the joint density function of ${\boldsymbol{X}}$ and $Z$ is

$f_{{\boldsymbol{X}}, Z}({\boldsymbol{x}}, z) = f_{{\boldsymbol{X}}|Z}({\boldsymbol{x}}, z)f_{Z}(z) = \frac{1}{ (2\pi)^{\frac{n}{2}}z^{\frac{n}{2}}|{\bf{\Sigma}}|^{\frac{1}{2}} } \exp \Big{\{} -\frac{1}{2z}{\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{x}}+{\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}} -\frac{z}{2}{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}-z 1_{ \{z\geq 0\} } \Big{\}}.$

Suppose that there are $T$ observations ${\boldsymbol{X}}_1, \ldots, {\boldsymbol{X}}_T$ generated from the latent random variables $z_1, z_2, \cdots, z_T$ respectively. The complete data is defined as $\{ ({\boldsymbol{x}}_t, z_t) \}\, , t = 1, 2, \cdots, T$. In the EM algorithm, ${\boldsymbol{x}}_t$ and $z_t$ are the observed and missing data respectively. The log-likelihood up to an additive constant can be written as

$\tilde{L}({\boldsymbol{\mu}}, {\bf{\Sigma}}) = \sum\limits_{t = 1}^{T}\ln f_{{\boldsymbol{X}}, Z}({\boldsymbol{x}}_t, z_t) = -\frac{T}{2}\ln|{\bf{\Sigma}}|-\frac{1}{2}\sum\limits_{t = 1}^{T}\frac{1}{z_t}{\boldsymbol{x}}'_t{\bf{\Sigma}}^{-1}{\boldsymbol{x}}_t +\sum\limits_{t = 1}^{T}{\boldsymbol{x}}'_t{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}} -\frac{1}{2}{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}\sum\limits_{t = 1}^{T}z_{t}-\sum\limits_{t = 1}^{T}z_{t}1_{ \{z_t\geq 0\} }.$

Note that the last term of the above equation does not contain any unknown parameters and thus is negligible. Then, the E-step becomes

$E\Big(\tilde{L}({\boldsymbol{\mu}}, {\bf{\Sigma}})|{\boldsymbol{x}}_{t}, \hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}}\Big) \propto -\frac{T}{2}\ln\big|{\bf{\Sigma}}\big| + \sum\limits_{t = 1}^{T}{\boldsymbol{x}}'_t{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}} - \frac{1}{2}\sum\limits_{t = 1}^{T}E(z_t^{-1}|{\boldsymbol{x}}_t, \hat{{\boldsymbol{\mu}}}, \hat\Sigma){\boldsymbol{x}}'_t{\bf{\Sigma}}^{-1}{\boldsymbol{x}}_t -\frac{1}{2}{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}\sum\limits_{t = 1}^{T}E(z_{t}|{\boldsymbol{x}}_t, \hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}}).$

where $E(z_t|{\boldsymbol{x}}, \hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}})$ and $E(z_t^{-1}|{\boldsymbol{x}}, \hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}})$ are the conditional expectations of $z_t$ and $z_t^{-1}$ given ${\boldsymbol{x}}_t$ and the current estimates $\hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}}$.

To evaluate conditional expectations $E(z_t^{-1}|{\boldsymbol{x}}_t, \hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}})$ and $E(z_t|{\boldsymbol{x}}_t, \hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}})$, we need the conditional density of $Z$ given ${\boldsymbol{X}}$, $f_{Z|{\boldsymbol{X}}}$. After some straightforward algebra, the conditional distribution of $Z$ given ${\boldsymbol{X}}$ is an inverse Gaussian distribution with density function

$\begin{eqnarray} f_{Z|{\boldsymbol{X}}}(z|{\boldsymbol{x}}, {\boldsymbol{\mu}}, {\bf{\Sigma}}) & = & \frac{f_{{\boldsymbol{X}}, Z}({\boldsymbol{x}}, z)}{f_{{\boldsymbol{X}}}({\boldsymbol{x}})} = \frac{ \frac{1}{ (2\pi)^{\frac{n}{2}}z^{\frac{n}{2}}|{\bf{\Sigma}}|^{\frac{1}{2}} } \exp \big{\{} -\frac{1}{2z}{\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{x}}+{\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}} -\frac{z}{2}{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}-z 1_{ \{z\geq 0\} } \big{\}} }{ \frac{ 2e^{{\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}} } { (2\pi)^{n/2} |{\bf{\Sigma}}|^{1/2} } \big{(} \frac{{\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{x}} }{ 2+{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}} \big{)} ^{v/2} K_{v} \big{(} \sqrt{ (2+{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}})({\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{x}}) } \big{)} } \nonumber \\ & = & \Big{(} \frac{ 2+{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}}{{\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{x}} } \Big{)} ^{v/2} z^{-n/2} \frac{ \exp \big{\{} -\frac{1}{2} \big{[} z^{-1} {\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{x}} +z({\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}+z 1_{\{z\geq 0\}}) \big{]} \big{\}} }{ 2K_{v} \big{(} \sqrt{ (2+{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}})({\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{x}}) } \big{)} }. \label{conditional-density} \end{eqnarray}$

(3.1)

Lemma 3.1. (GIG (Stacy, 1962)) A random variable $X$ follows Generalized Inverse Gaussian distribution(denoted as $X\sim N^{-}(\lambda, \chi, \psi)$) if its density function could be represented as

$f(x) = \frac{ \chi^{-\lambda}(\sqrt{\chi\psi})^{\lambda} }{2K_{\lambda}(\sqrt{\chi\psi})} x^{\lambda-1} \exp\Big{\{}-\frac{1}{2}(\chi x^{-1}+\psi x)\Big{\}} \, , \quad x>0.$

where $K_{\lambda}$ denotes the third kind modified Bessel function, and the parameters satisfy

$\begin{equation*} \left\{ \begin{array}{l} \chi>0, \psi\geq 0, & if \;\lambda<0; \\ \chi>0, \psi>0, & if \;\lambda = 0; \\ \chi\geq0, \psi>0, & if\; \lambda>0. \end{array} \right. \end{equation*}$

After some algebraic manipulations, it is easy to show that $Z|{\boldsymbol{X}}$ follows Generalized Inverse Gaussian distribution:

$Z|{\boldsymbol{X}} \sim N^{-}\Big{(} \frac{2-n}{2}, {\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{x}}\, , 2+{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}} \Big{)}.$

If $\chi>0, \psi>0,$ the moments could be calculated through the following formulas:

$\begin{eqnarray*} && E(X^{\alpha}) = \Big{(} \frac{\chi}{\psi} \Big{)}^{\alpha/2} \frac{K_{\lambda+\alpha}(\sqrt{\chi\psi})}{K_{\lambda}(\sqrt{\chi\psi})}, \quad \alpha\in \mathbb{R}, \\ && E(\ln X) = \frac{dE(X^{\alpha})}{d\alpha}\Big{|}_{\alpha = 0}. \end{eqnarray*}$

Denote $\chi = {\boldsymbol{x}}'{\bf{\Sigma}}^{-1}{\boldsymbol{x}}, \psi = 2+{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}$. Then, $Z|{\boldsymbol{X}}\sim N^-(\frac{2-n}{2}, \chi, \psi)$. From the conditional density function of (3.1), we can obtain the conditional expectations with the following moment properties:

$\begin{eqnarray*} a_{t} & = & E(z_t|{\boldsymbol{x}}_t, \hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}}) = \sqrt{\frac{\chi_t}{\psi}}\frac{K_{\frac{n}{2}-2}(\sqrt{\chi_t\psi})}{K_{\frac{n}{2}-1}(\sqrt{\chi_t\psi})}\, , \quad t = 1, 2, \cdots, T;\\ b_{t} & = & E(z_{t}^{-1}|{\boldsymbol{x}}_t, \hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}}) = \sqrt{\frac{\psi}{\chi_t}}\frac{K_{\frac{n}{2}}(\sqrt{\chi_t\psi})}{K_{\frac{n}{2}-1}(\sqrt{\chi_t\psi})}\, , \quad t = 1, 2, \cdots, T. \end{eqnarray*}$

where $\chi_t = {\boldsymbol{x}}'_t\Sigma^{-1}{\boldsymbol{x}}_t$, $R(\lambda) = \frac{ K_{\lambda+1}(x) }{ K_{\lambda}(x) }$ is strictly decreasing in $x$ with $\lim_{x\rightarrow \infty}R_{\lambda}(x) = 1$ and $\lim_{x\rightarrow 0^+}R_{\lambda}(x) = \infty$. Thus, $a_t>0, b_t>0\, , t = 1, 2, \cdots, T$.

Finally, if the conditional expectation $E(z_t|{\boldsymbol{x}}_t, \hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}})$ and $E(z_t^{-1}|{\boldsymbol{x}}_t, \hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}})$ are replaced by $a_t$ and $b_t$ respectively, the objective function becomes

$Q({\boldsymbol{\mu}}, {\bf{\Sigma}}|{\boldsymbol{x}}_t, \hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}}) = -\frac{T}{2}\ln|{\bf{\Sigma}}|+\sum\limits_{t = 1}^{T}{\boldsymbol{x}}'_t{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}} -\frac{1}{2}\sum\limits_{t = 1}^{T}b_{t}{\boldsymbol{x}}'_t{\bf{\Sigma}}^{-1}{\boldsymbol{x}}_t -\frac{1}{2}{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}\sum\limits_{t = 1}^{T}a_{t}.$

(3.2)

Denote ${\boldsymbol{S}} = {\bf{\Sigma}}^{-1}$. The objective function (3.2) becomes

$Q({\boldsymbol{\mu}}, {\bf{\Sigma}}|{\boldsymbol{x}}_t, \hat{{\boldsymbol{\mu}}}, \hat{{\boldsymbol{S}}}) = \frac{T}{2}\ln|{\boldsymbol{S}}|+\sum\limits_{t = 1}^{T}{\boldsymbol{x}}'_t{\boldsymbol{S}}{\boldsymbol{\mu}} -\frac{1}{2}\sum\limits_{t = 1}^{T}b_{t}{\boldsymbol{x}}'_t{\boldsymbol{S}}{\boldsymbol{x}}_t -\frac{1}{2}{\boldsymbol{\mu}}'{\boldsymbol{S}}{\boldsymbol{\mu}}\sum\limits_{t = 1}^{T}a_{t}.$

(3.3)

Taking derivative of objective function (3.3) with respect to ${\boldsymbol{\mu}}, {\boldsymbol{S}}$, we obtain

$\begin{eqnarray*} \frac{\partial Q}{\partial {\boldsymbol{\mu}}} & = & \sum\limits_{t = 1}^{T} {\boldsymbol{x}}'_t {\boldsymbol{S}} - \sum\limits_{t = 1}^{T} a_t {\boldsymbol{\mu}}'{\boldsymbol{S}} = 0, \\ \frac{\partial Q}{\partial {\boldsymbol{S}}} & = & \frac{T}{2}{\boldsymbol{S}}^{-1}-\frac{1}{2}\sum\limits_{t = 1}^{T} b_t {\boldsymbol{x}}'_t {\boldsymbol{x}}_t+\sum\limits_{t = 1}^{T}{\boldsymbol{x}}_t'{\boldsymbol{\mu}} -\frac{1}{2}\sum\limits_{t = 1}^{T} a_t{\boldsymbol{\mu}}'{\boldsymbol{\mu}} = 0. \end{eqnarray*}$

Substituting ${\boldsymbol{S}}$ by ${\bf{\Sigma}}$ and setting these derivatives to zero yield

$\begin{eqnarray*} && \sum\limits_{t = 1}^{T} {\boldsymbol{x}}'_t{\bf{\Sigma}}^{-1}- \sum\limits_{t = 1}^{T} a_t {\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1} = {\boldsymbol{0}}, \\ && \frac{T}{2}{\bf{\Sigma}}-\frac{1}{2}\sum\limits_{t = 1}^{T}b_{t}{\boldsymbol{x}}'_t {\boldsymbol{x}}_t+\sum\limits_{t = 1}^{T}{\boldsymbol{x}}'_t{\boldsymbol{\mu}} -\frac{1}{2}\sum\limits_{t = 1}^{T}a_{t}{\boldsymbol{\mu}}'{\boldsymbol{\mu}} = {\boldsymbol{0}}. \end{eqnarray*}$

Thus, maximization of the objective function $Q({\boldsymbol{\mu}}, {\bf{\Sigma}}|{\boldsymbol{x}}_t, \hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}})$ can be achieved by the following iterative updating formulas:

$\hat{{\boldsymbol{\mu}}} = \frac{\bar{{\boldsymbol{x}}}}{\bar{a}}; \quad \hat{{\bf{\Sigma}}} = \overline{b_{t}{\boldsymbol{x}}'_t {\boldsymbol{x}}_t}-\frac{ \bar{{\boldsymbol{x}}}'\bar{{\boldsymbol{x}}} }{\bar{a}}.$

where $\bar{a}, \bar{b}$ stand for the average of $\{a_t\}_{t = 1}^T$ and $\{b_t\}_{t = 1}^T$ respectively and $\bar{{\boldsymbol{x}}}$ is the average of $\{{\boldsymbol{x}}_t\}_{t = 1}^{T}$. In what follows, we present the iterative reweighted Expectation-Maximization algorithm for parameter estimation of Asymmetric Laplace distribution.

Algorithm 1 Iterative reweighting algorithm

1. Set iteration number $k=1$ and select initial estimates for parameters ${\boldsymbol{\mu}}^{(0)}, {\bf{\Sigma}}^{(0)}$. 2. (E-Step) At $k$-th iteration with current estimates ${\boldsymbol{\mu}}^{(k)}, {\bf{\Sigma}}^{(k)}$, define the corresponding log-likelihood as                         $l^{(k)} = \log \sum\limits_{t=1}^{T} f({\boldsymbol{x}}_t|{\boldsymbol{\mu}}^{(k)}, {\bf{\Sigma}}^{(k)}), \quad k=1, 2, \cdots.$ With notations $\chi_t ={\boldsymbol{x}}'_t{\bf{\Sigma}}^{-1}{\boldsymbol{x}}_t, \psi = 2+{\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}$, we can obtain iterative weights                         $a_{t} = E(z_t|{\boldsymbol{x}}_t, \hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}}) =\sqrt{\frac{\chi_t}{\psi}}\frac{K_{2-\frac{n}{2}}(\sqrt{\chi_t\psi})}{K_{1-\frac{n}{2}}(\sqrt{\chi_t\psi})}, \quad t=1, 2, \cdots, T;$                         $b_{t} = E(z_{t}^{-1}|{\boldsymbol{x}}_t, \hat{{\boldsymbol{\mu}}}, \hat{{\bf{\Sigma}}}) =\sqrt{\frac{\psi}{\chi_t}}\frac{K_{-\frac{n}{2}}(\sqrt{\chi_t\psi})}{K_{1-\frac{n}{2}}(\sqrt{\chi_t\psi})}, \quad t=1, 2, \cdots, T.$ 3. (M-Step) Employ the following iteration formulas to calculate the new estimates ${\boldsymbol{\mu}}^{(k+1)}, {\bf{\Sigma}}^{(k+1)}$ at $(k+1)$-th iteration:                         ${\boldsymbol{\mu}}^{(k+1)}=\frac{\bar{{\boldsymbol{x}}}}{\bar{a}}, \quad {\bf{\Sigma}}^{(k+1)} = \overline{b_{t}{\boldsymbol{x}}'_t {\boldsymbol{x}}_t}-\frac{ \bar{{\boldsymbol{x}}}'\bar{{\boldsymbol{x}}} }{\bar{a}}.\quad \quad \quad \quad (3.4)$ The log-likelihood at $(k+1)$-th iteration becomes                         $l^{(k+1)} = \log \sum\limits_{t=1}^{T} f({\boldsymbol{x}}_t|{\boldsymbol{\mu}}^{(k+1)}, {\bf{\Sigma}}^{(k+1)}).$ 4. Repeat these iteration steps until convergence with criterion $l^{(k+1)}-l^{(k)} < \varepsilon$, where $\varepsilon>0$ is a small number that control the convergence precision, for convenience, we take $\varepsilon=1e^{-16}$.

---

4. Simulation studies

To evaluate the performance of portfolio selection models and parameter estimation methods in Section 3, we generate 100 datasets from Gaussian distribution and Asymmetric Laplace distribution respectively. Each dataset consists of $T = 200$ observations with the following parameter settings: Case (1): $n = 3\, , {\boldsymbol{\mu}} = (0.03\, , 0.06\, , 0.09)$; Case (2): $n = 5\, , {\boldsymbol{\mu}} = (0.01\, , 0.02\, , 0.06\, , 0.08\, , 0.09)$; Case (3): $n = 10\, , {\boldsymbol{\mu}} = (0.01\, , 0.02\, , 0.03\, , \cdots, 0.10)$. For each case, we set ${\bf{\Sigma}} = \text{diag}\, ({\boldsymbol{\mu}}/10)$. All the simulation studies are carried out on a PC with Intel Core i7 3.6 GHz processor under $\texttt{R}$ platform.

Each dataset are estimated under both multivariate Gaussian and Asymmetric Laplace distribution. ALD (EM-AL) is estimated using the EM algorithm described in Section 3.2. We evaluate the estimation performance using Bias measure, defined as $\text{Bias} = \|\hat{{\boldsymbol{\mu}}}-{\boldsymbol{\mu}}\|_1 + \|\hat{{\bf{\Sigma}}}-{\bf{\Sigma}} \|_1$. The mean log-likelihood and mean bias of the simulated 200 datasets are reported in Table 1.

Table 1 indicate that if the model is correctly specified, the estimation performance is always the best in terms of bias. If the data is generated from Gaussian distribution, the estimation from Gaussian model is the best, so does Asymmetric Laplace distribution. If data is generated from Gaussian distribution, then the estimation log-likelihood of Gaussian model is larger than Asymmetric Laplace distribution, this is true for Asymmetric Laplace data as well.

Figure 3 show that for Gaussian data, since Gaussian data fit the model better, efficient frontiers under Gaussian data are more close to Gaussian models; Figure 4 indicate that for generated Asymmetric Laplace data, efficient frontiers nearly equivalent to true Asymmetric Laplace data. Figure 3–4 suggest that we can first modeling data using Gaussian and Asymmetric Laplace distribution, and use the fitted log-likelihood to determine the distribution, then we evaluate the performance with the corresponding efficient frontier analysis.

---

5. Real data analysis

In this section, we apply our proposed methodology to two real financial datasets, Hang Seng Index and Nasdaq Index, both datasets are downloaded from Bloomberg, with daily data range from January 4, 2011, to December 29, 2017. The variable of interest is the rate of returns multiplied by the annualized ratio $\sqrt{252}$, formulated as

$\text{LogRet}\, (t) = \sqrt{252} \Big\{ \log\big(\text{price}[t+1]\big)-\log\big(\text{price}[t]\big) \Big\}, \quad t = 1, 2, \cdots, 1721.$

These two datasets are analyzed through efficient frontier analysis under ALD framework on $\texttt{R}$ platform. Theorem 3.1–3.2 indicate that portfolio selection models under ALD framework can be reduced to the following quadratic programming problem:

$\min\limits_{{\boldsymbol{w}}} \sigma^2 = {\boldsymbol{w}}' {\bf{\Sigma}} {\boldsymbol{w}} \quad \text{ s.t. } \quad {\boldsymbol{w}}'{\boldsymbol{\mu}} = r_0\, , {\boldsymbol{w}}'{\boldsymbol{1}} = 1.$

with explicit solution (see Lai and Xing, 2008) as follows:

$\hat{{\boldsymbol{w}}} = \frac{D-r_0B}{AD-B^2}{\bf{\Sigma}}^{-1}{\boldsymbol{1}} + \frac{r_0A-B}{AD-B^2}{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}.$

(5.1)

Here, $A = {\bf{1}}'{\bf{\Sigma}}^{-1}{\bf{1}}, B = {\bf{1}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}$ and $D = {\boldsymbol{\mu}}'{\bf{\Sigma}}^{-1}{\boldsymbol{\mu}}$.

---

5.1. Example 1: Hang seng index

In the first example, we construct a portfolio consisting of 8 Hang Seng indexes: HK1, HK175, HK2007, HK2318, HK4, HK6, HK66. The summary descriptive statistics are reported in Table 2. It is clear that the all stock returns exhibit larger skewness and kurtosis. The median of these stocks are close to zero, the log-likelihood of Asymmetric Laplace distribution is larger than gaussian distribution, indicating that Asymmetric Laplace distribution would be a good fit than gaussian distribution.

Then we fit the data to Asymmetric Laplace distribution through EM algorithm as described in Section 3.2. Parameter estimation results are displayed in Table 2, we construct portfolios under Asymmetric Laplace framework at different levels of expected return. Consider increasing target expected return values

$r_0 = 0.040\, , 0.0745\, , 0.1081\, , 0.1417\, , 0.1753\, , 0.2088\, , 0.2424\, , 0.2760\, , 0.3096\, , 0.3431\, .$

Portfolio selection results are summarized in Table 3, with kurtosis, skewness, sharpe ratio, VaR and ES results at $\alpha = 0.01, 0.05, 0.10$ levels.

The efficient frontier tendencies are displayed in Figure 5. It is suggested that aggressive investors should impose higher confidence levels and conservative investors may choose smaller confidence levels. Figure 6 depicts the kurtosis, skewness, and sharpe ratio tendency of portfolio selection models. Results show that Sharpe Ratio, Skewness and Kurtosis increase fast and decreases slowly down as the target expected returns increases.

---

5.2. Example 2: Nasdaq index

In the second example, we consider Nasdaq index, including CTRP, MNST, NFLX, NTES, NVDA, TTWO, and report the descriptive statistics in Table 4. All the indexes exhibit significant skewness and kurtosis. Jarq.Test results indicate that this dataset deviates from normality significantly. We fit the log returns data to Gaussian and Asymmetric Laplace distributions. Since Asymmetric Laplace model achieve higher log-likelihood results compared to Gaussian model, we choose EM-AL model for data fitting. Parameter estimation results are displayed in Table 4.

Then we consider increasing target expected returns

$r_0 = 0.2100, 0.4958, 0.7815, 1.0673, 1.3530, 1.6388, 1.9245, 2.2102, 2.4960, 2.7817.$

Results of skewness, kurtosis, sharpe ratio and VaR, ES results are summarized in Table 5. Figure 7 displays the efficient frontiers at confidence level $\alpha = 0.01, 0.05, 0.10$. These results show that the portfolio capture higher risk at higher $\alpha$ levels. Figure 8 displays the skewness, kurtosis, and Sharpe Ratio Tendency. The optimal portfolios can be obtained from (5.1) with the corresponding VaR, ES, skewness, kurtosis and Sharpe Ratio.

Table 5 and Figure 7 suggests that as $r_0$ increases, all ES (ES$_{0.01}$, ES$_{0.05}$, ES$_{0.10}$) increases, indicating that higher return is derived from higher risk. It is interesting that under the ALD assumption, as $r_0$ increases, Sharpe ratio and skewness first decreases then increases accordingly. As $\alpha$ increases, VaR and ES measures decreases. Thus, conservative investors can choose larger $\alpha$ levels and aggressive investors would select smaller $\alpha$ levels.

---

6. Conclusion and prospects

In this paper, we derive several equivalent portfolio selection models under ALD framework, these models can be transformed to quadratic programming problem with explicit solutions. The Expectation-Maximization algorithm for parameter estimation of Asymmetric Laplace distribution is obtained and outperforms moment estimation.

There are several advantages of Asymmetric Laplace distribution based models. First, the equivalence of risk measures such as VaR, ES and StD at maximization of skewness and minimization of kurtosis faciliate portfolio selection models significantly. Second, confidence levels of these models offer investors various portfolio selection choices. Conservative investors can choose larger $\alpha$ levels and aggressive investors can select smaller $\alpha$ levels. Finally, ALD model is able to explain skewness and kurtosis in financial data. Therefore, the Asymmetric Laplace distribution can be widely applied to handle real financial datasets.

We may further extend the Asymmetric Laplace based portfolio selection model to cases of mixture Asymmetric Laplace distributions. Another direction is to combine clustering techniques (see Dias et al., 2015; Iorio et al., 2018) with Asymmetric Laplace distribution for portfolio selection of time series models.

---

Acknowledgments

Chi Tim, Ng's work is supported by the 2016 Chonnam National University Research Program grant (No. 2016–2762).

---

Conflict of interest

The authors declare no conflict of interest.

---