#  Non-parametric cumulants approach for outlier detection of multivariate
financial data

Francesco Cesarone1, Rosella Giacometti2, Jacopo Maria Ricci3
1 _Roma Tre University - Department of Business Studies_
<EMAIL_ADDRESS>
2 _University of Bergamo - Department of Management_
<EMAIL_ADDRESS>
2 _University of Bergamo - Department of Economics_
<EMAIL_ADDRESS>

###### Abstract

In this paper, we propose an outlier detection algorithm for multivariate data
based on their projections on the directions that maximize the Cumulant
Generating Function (CGF). We prove that CGF is a convex function, and we
characterize the CGF maximization problem on the unit $n$-circle as a concave
minimization problem. Then, we show that the CGF maximization approach can be
interpreted as an extension of the standard principal component technique.
Therefore, for validation and testing, we provide a thorough comparison of our
methodology with two other projection-based approaches both on artificial and
real-world financial data. Finally, we apply our method as an early detector
for financial crises.

Keywords: Outlier detection; Cumulant generating function; Principal component
analysis; Projections; Multivariate skew-normal distribution.

##  1 Introduction

Outlier detection issue has become increasingly important over the years and,
as of now, its fields of application range from medicine and engineering to
finance (see, e.g., Singh and Upadhyaya, 2012; Meng et al., 2019). As for the
latter, outliers can be the consequence of human error or fraudulent
activities; similarly, financial crises can be viewed as anomalies since
markets experience atypical behaviors in those periods (see, e.g., Ané et al.,
2008). Furthermore, a few outliers can strongly influence the results of an
experiment. This can be observed in Portfolio Optimization, where some
portfolio selection models can be highly sensitive to changes in input data
(see, e.g., Kondor et al., 2007; Cesarone et al., 2020; Giacometti et al.,
2021). Because of this widespread practical relevance, many authors tackled
this topic. Hence, the theory behind anomaly detection has unsurprisingly
evolved, from the first studies which dealt with more simple instances, i.e.,
univariate Gaussian data, to more complex cases, such as multivariate data
following non-parametric distributions. We can mention the work of Ferguson
(1961), who considers univariate normal samples and identifies outliers as
data with mean slippage. Wilks (1963) proposes a method to identify outliers
in a multivariate normal distribution with unknown parameters. Gnanadesikan
and Kettenring (1972) propose outlier detection methods on multivariate data,
based on their projection onto the directions corresponding to the principal
components obtained by the standard Principal Component Analysis (PCA).
Schwager and Margolin (1982) extend the work of Ferguson (1961) to
multivariate normal data. Reed and Yu (1990) identify outliers by using the
so-called RX detector. Such an indicator measures the location of multivariate
data points in the dispersion ellipsoid by means of the Mahalanobis distance.
The authors assume that both ordinary and outlier data follow a multivariate
normal distribution with same covariance matrix but different mean vectors.
Note that, even though Reed and Yu (1990) assume the Gaussian distribution,
other authors used the RX detector, removing such a strong distributional
hypothesis. For example, Das and Sinha (1986) and Sinha (1984) extend the work
by Schwager and Margolin (1982) to (nonnormal) elliptically symmetric
distributions, and use the Mahalanobis distance as a detector. Since outliers
affect both location and scale parameters, one way to overcome this issue is
to use robust estimates. For example, Maronna (1976) proposes robust
M-estimators of multivariate mean vectors and covariance matrices. Both Stahel
(1981) and Donoho (1982) define, in their respective works, an affine
equivariant robust estimator for location and scatter. Another popular method
consists in computing the minimal covariance determinant (see, e.g.,
Rousseeuw, 1984, 1985).

Especially when dealing with high-dimensional multivariate data, many
techniques aim to find outliers in the univariate projections of such data to
reduce the computational effort. For this reason, several studies have been
devoted to identifying the “best” directions in which data must be projected
to represent the variability of data most effectively. In Peña and Prieto
(2001, 2007), these directions correspond to those that maximize and minimize
the kurtosis of the projected data. Domino (2020) generalizes the Prieto and
Peña’s approach by choosing the direction that maximizes the fourth cumulant.
Following this stream of literature, our work aims at detecting outliers, by
projecting the data onto the direction that maximizes the Cumulant Generating
Function (CGF).

In this paper, we refine some theoretical results of the methodologies
proposed by Bernacchia and Naveau (2008) and Bernacchia et al. (2008). More
precisely, we prove that CGF is a convex function, and, then, we characterize
the CGF maximization problem on the unit $n$-circle as a concave minimization
problem. Then, we extend the outlier detection methodology, based on the
projections of multivariate data on the directions obtained by the classical
PCA technique, to the directions that maximize CGF. Finally, we perform an
extensive empirical analysis both on simulated and historical data, comparing
our method with those described in Peña and Prieto (2001, 2007) and Domino
(2020).

The paper is organized as follows. In Section 2, we first introduce some
preliminary concepts about the moment generating function, the cumulant
generating function, and cumulants for univariate and multivariate random
variables. Then, we report and refine some results developed in Bernacchia and
Naveau (2008) and Bernacchia et al. (2008) for a generalization of the
principal component analysis (PCA) technique. In Section 3, we describe our
outlier detection algorithm, and the methods used to comparative purposes. For
validation and testing, in Section 4, we present a thorough comparison of
these methodologies both on artificial and real-world data. The real-world
data application identifies the outliers as the materialization of financial
crises. Finally, Section 5 contains some concluding remarks.

##  2 Theoretical framework

For the sake of completeness and readability, we recall below some notions
about the moment generating function, the cumulant generating function, and
the cumulants in the case both of a univariate and a multivariate random
variable. Furthermore, we also report and refine some concepts introduced by
Bernacchia and Naveau (2008) and Bernacchia et al. (2008) for a generalization
of the principal component analysis (PCA) technique.

Let $X$ be a univariate random variable. If $\mathbb{E}[e^{\xi X}]$ exists and
is finite $\forall\xi\in\mathbb{R}$, then the moment generating function of
$X$ is defined as follows

$M_{X}(\xi)=\mathbb{E}[e^{\xi
X}]=\sum_{m=0}^{+\infty}\mu_{m}\frac{\xi^{m}}{m!}$ (1)

where $\mu_{m}=\mathbb{E}[X^{m}]$ is the $m^{th}$ raw moment of $X$, that can
be obtained by differentiating $m$ times w.r.t. $\xi$ and setting $\xi=0$,
namely $\mu_{m}=M_{X}^{(m)}(0)$. The cumulant generating function (CGF) of $X$
can be expressed as the logarithm of (1)

$G_{X}(\xi)=\ln\mathbb{E}[e^{\xi
X}]=\sum_{m=1}^{+\infty}k_{m}\frac{\xi^{m}}{m!},$ (2)

where $k_{m}$ is the $m^{th}$ order cumulant of $X$. It is straightforward to
see that $k_{m}$ can be obtained by differentiating Expression (2) $m$ times
and setting $\xi=0$, namely $k_{m}=G_{X}^{(m)}(0)$. Furthermore, the cumulants
can be expressed as functions of the moments. For instance, the first 4
cumulants of $X$ are as follows

$\begin{array}[]{l}k_{1}=\mathbb{E}[X]\\\
k_{2}=\mathbb{E}[(X-\mathbb{E}[X])^{2}]\\\
k_{3}=\mathbb{E}[(X-\mathbb{E}[X])^{3}]\\\
k_{4}=\mathbb{E}[(X-\mathbb{E}[X])^{4}]-3(\mathbb{E}[(X-\mathbb{E}[X])^{2}])^{2}\\\
\end{array}$

Note that the first two cumulants correspond to the mean and variance
respectively, the third coincides with the third central moment, while from
the fourth onward the cumulants are polynomial functions of the central
moments with integer coefficients. From Expression (2), we also observe that
for small values of $\xi$ the cumulant generating function of $X$ is
essentially determined by its variance (if $k_{1}=0$), whereas for large $\xi$
the contributions of the higher order cumulants become dominant. Furthermore,
if $X\sim N(\mu,\sigma^{2})$, then
$G_{X}(\xi)=\mu\xi+\sigma^{2}\displaystyle\frac{\xi^{2}}{2}$ since $k_{m}=0$
for $m>2$.

In the case of a multivariate random variable
$\boldsymbol{X}=(X_{1},\ldots,X_{n})$, its CGF becomes

$G_{\boldsymbol{X}}(\boldsymbol{\xi})=\ln\mathbb{E}[e^{\boldsymbol{\xi}^{T}\boldsymbol{X}}],$
(3)

where $\boldsymbol{\xi}=(\xi_{1},\ldots,\xi_{n})\in\mathbb{R}^{n}$. Denoting
by $r=\|\boldsymbol{\xi}\|_{2}$ the euclidian norm of $\boldsymbol{\xi}$ and
by $\boldsymbol{\theta}$ the versor of $\boldsymbol{\xi}$, we can write
$\boldsymbol{\xi}=r\boldsymbol{\theta}$, and therefore the cumulant generating
function of $\boldsymbol{X}$ can be expressed as follows

$G_{\boldsymbol{X}}(r,\boldsymbol{\theta})=\ln\mathbb{E}[e^{r\boldsymbol{\theta}^{T}\boldsymbol{X}}]=\sum_{m=1}^{+\infty}k_{m}(\boldsymbol{\theta})\frac{r^{m}}{m!}.$
(4)

where
$k_{m}(\boldsymbol{\theta})=\displaystyle\frac{d^{m}G_{\boldsymbol{X}}(r,\boldsymbol{\theta})}{dr^{m}}\biggl{\lvert}_{r=0}$
is the $m^{th}$ cumulant of $\boldsymbol{X}$ projected along the direction
$\boldsymbol{\theta}$.

Note that the length $r$ of $\boldsymbol{\xi}$ plays the same rule of $\xi$ in
the univariate case, namely if $r$ is small then the information on
$\boldsymbol{X}$ through CGF is basically represented by the covariance matrix
of $\boldsymbol{X}$ (if data are centered, i.e., $\mathbb{E}[X_{j}]=0$ with
$j=1,\ldots,n$). Whereas if $r$ is large then the information about
$\boldsymbol{X}$ described by its CGF predominantly depends on the higher-
order cumulants.

###  2.1 PCA through CGF

Following the work of Bernacchia and Naveau (2008), we look for identifying
the largest principal components, namely the directions that provide most of
the variability present in the original data, by exploiting the information
contained in the Cumulant Generating Function (CGF).

Using the classical PCA technique the principal component with the largest
eigenvalue is the versor maximizing the variance, which, as we will show in
the next section, also coincides with the optimal versor maximizing CGF (4) in
the case of multivariate normal random variables, or in the case of generic
random vectors but for small values of $r$. As mentioned by Bernacchia and
Naveau (2008), the CGF maximization aims to find the directions with the
largest variability of the multivariate distribution not only through the
first two cumulants, but also through the higher-order ones, which become
dominant especially for data points that deviate widely from the mean of the
random vector. Hence, for a fixed $r$, our goal is to find the directions that
maximize (4), i.e.

$\begin{array}[]{lll}\displaystyle\max_{\boldsymbol{\theta}}&G_{\boldsymbol{X}}(r,\boldsymbol{\theta})\\\
\mbox{s.t.}&&\\\ &\boldsymbol{\theta}^{T}\boldsymbol{\theta}=1&\end{array}$
(5)

As we will show in Section 2.3, Problem (5) is a constrained concave
programming problem, where the main difficulty lies in the fact that concave
problems normally have many local maxima points. In the next section, we
present our approach on two artificial multivariate random vectors, normal and
skew-normal.

###  2.2 Some special cases: multivariate normal and skew-normal

To better understand the rationale behind the approach of Bernacchia and
Naveau (2008), we report the CGF maximization procedure in the case of
multivariate normal (Section 2.2.1) and skew-normal (Section 2.2.2) random
variables, while in Section 2.3 we address the general case. We show that, in
the case of a multivariate normal random vector, the optimal versor maximizing
CGF collapses to the first principal component of the PCA. For a skew-normal,
we will see that the optimal solution obtained from Problem (5) tends to the
first principal component when the radius of the hypersphere is small, while,
in the case of a large radius, the direction of the optimal versor deviates
from the first component toward the direction where the multivariate
distribution presents a fatter tail.

Without loss of generality, we henceforth assume that the data are centered
around their mean.

####  2.2.1 Multivariate normal case

Let $\boldsymbol{X}\sim N(\boldsymbol{0},\boldsymbol{\Sigma})$, where,
therefore, $\boldsymbol{X}$ is mean-centered and $\boldsymbol{\Sigma}$ is its
covariance matrix. Thus, the cumulant generating function of $\boldsymbol{X}$
is

$G_{\boldsymbol{X}}(\xi)=\frac{1}{2}\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}\,.$

Indeed, similar to the univariate case, for normal random vectors, the first
two cumulants are exactly the mean and variance, while the higher-order
cumulants are all equal to zero. Then, Problem (5) can be rewritten as

$\begin{array}[]{lll}\displaystyle\max_{\boldsymbol{\theta}}&\displaystyle\frac{r^{2}}{2}\boldsymbol{\theta}^{T}\boldsymbol{\Sigma}\boldsymbol{\theta}&\\\
\mbox{s.t.}&&\\\ &\boldsymbol{\theta}^{T}\boldsymbol{\theta}=1&\end{array}$
(6)

This problem can be easily solved by means of the Lagrange multipliers method,
where the Lagrangian function is
$L(\lambda,\boldsymbol{\theta})=\displaystyle\frac{r^{2}}{2}\boldsymbol{\theta}^{T}\boldsymbol{\Sigma}\boldsymbol{\theta}-\lambda(\boldsymbol{\theta}^{T}\boldsymbol{\theta}-1)$.
Hence, the corresponding first order conditions are
$\nabla_{\lambda}L(\lambda,\boldsymbol{\theta})=\boldsymbol{\theta}^{T}\boldsymbol{\theta}-1=0$
and
$\nabla_{\boldsymbol{\theta}}L(\lambda,\boldsymbol{\theta})=r^{2}\boldsymbol{\Sigma}\boldsymbol{\theta}-2\lambda\boldsymbol{\theta}=0$.
Problem (6) has two optimal solutions, the versors
$\boldsymbol{\hat{\theta}_{1}}$ and
$\boldsymbol{\hat{\theta}_{2}}=-\boldsymbol{\hat{\theta}_{1}}$, which are both
eigenvectors of $\boldsymbol{\Sigma}$. Note that applying the classical PCA
technique, the first principal components can be obtained by solving the
following problem

$\begin{array}[]{lll}\displaystyle\max_{\boldsymbol{\theta}}&\boldsymbol{\theta}^{T}\boldsymbol{\Sigma}\boldsymbol{\theta}\\\
\mbox{s.t.}&&\\\ &\boldsymbol{\theta}^{T}\boldsymbol{\theta}=1&\end{array}$
(7)

Since Problems (6) and (7) are equivalent except for the parameter
$\displaystyle\frac{r^{2}}{2}$ in the objective function of (6), these
problems lead to the same optimal solutions.

Figure 1: Bivariate normal random variable: the dashed lines represent its
isoprobability contours of $\boldsymbol{X}$, the red solid lines are the
isocurves of
$G_{\boldsymbol{X}}(r,\boldsymbol{\theta})=\frac{1}{2}\boldsymbol{\xi}^{\prime}\boldsymbol{\Sigma}\boldsymbol{\xi}$,
and the gray solid line is the circle of radius 1.

In Fig. 1, we give a graphical representation of the CGF maximization
procedure in the case of a bivariate normal random variable. The red solid
lines are the isocurves of
$G_{\boldsymbol{X}}(\boldsymbol{\xi})=\frac{1}{2}\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}$,
the gray solid line is the ball of radius 1, and the dashed lines represent
the isoprobability contours of $\boldsymbol{X}\sim
N(\boldsymbol{0},\boldsymbol{\Sigma})$. The points, where the (red solid)
isocurves and the (gray solid) cicle of radius 1 are tangent, represent the
optimal solutions of Problem (6), which coincide (except for a minus sign)
with the first principal component lying on the major axis of the (dashed)
local-dispersion ellipsoids, i.e., the direction of maximum variability of the
bivariate random variable. We point out that the fact that
$\boldsymbol{\hat{\theta}_{1}}$ and $\boldsymbol{\hat{\theta}_{2}}$ lie on the
same direction is due to the symmetry of the bivariate normal distribution.
Indeed, as we will see in the next section, in the case of a non-symmetric
distribution $\boldsymbol{\hat{\theta}_{1}}$ and
$\boldsymbol{\hat{\theta}_{2}}$ are in general on different directions, and
these directions depend on $r$. Clearly, this phenomenon will be more evident
for large $r$ and for highly skewed distributions.

####  2.2.2 Multivariate skew-normal case

We report here the implementation of the CGF maximization procedure in the
case of a multivariate skew-normal random variable.

Let $\boldsymbol{X}\sim
SN(\boldsymbol{\eta},\boldsymbol{\Sigma},\boldsymbol{\alpha})$, where
$\boldsymbol{\eta}$, $\boldsymbol{\Sigma}$ and $\boldsymbol{\alpha}$ are the
location, the scale and the shape parameters, respectively. As discussed in
Azzalini and Valle (1996), the probability density function (pdf) of
$\boldsymbol{X}$ is

$f_{\boldsymbol{X}}(\boldsymbol{x})=2\phi_{n}(\boldsymbol{x}-\boldsymbol{\eta};\boldsymbol{\Sigma})\Phi(\boldsymbol{\alpha}^{T}\operatorname{diag}(\boldsymbol{\sigma})^{-1}\bigl{(}\boldsymbol{x}-\boldsymbol{\eta}\bigr{)})\quad\mbox{with}\quad\boldsymbol{x}\in\mathbb{R}^{n},$
(8)

where $\phi_{n}(\boldsymbol{x}-\boldsymbol{\eta};\boldsymbol{\Sigma})$ is the
pdf of an $n$-variate Gaussian random variable with mean $\boldsymbol{\eta}$
and covariance matrix $\boldsymbol{\Sigma}$, $\Phi(\cdot)$ is the cumulative
distribution function (cdf) of a univariate standard normal random variable,
$\boldsymbol{\alpha}$ is the degree of skewness (when
$\boldsymbol{\alpha}=\boldsymbol{0}$, (8) collapses to the normal case), and
$\operatorname{diag}(\boldsymbol{\sigma})=\operatorname{diag}(\sigma_{1},\ldots,\sigma_{n})$,
i.e., it represents the diagonal matrix of the standard deviations. As shown
in Arellano-Valle and Azzalini (2008) and Azzalini and Capitanio (1999), we
have

$\boldsymbol{\mu_{X}}=\mathbb{E}[\boldsymbol{X}]=\boldsymbol{\eta}+\operatorname{diag}(\boldsymbol{\sigma})\boldsymbol{\delta}$
(9)
$\boldsymbol{\Sigma_{X}}=\mathbb{V}\mathrm{ar}[\boldsymbol{X}]=\boldsymbol{\Sigma}-\operatorname{diag}(\boldsymbol{\sigma})\boldsymbol{\delta}\boldsymbol{\delta}^{T}\operatorname{diag}(\boldsymbol{\sigma})$
(10)
$G_{\boldsymbol{X}}(\boldsymbol{\xi})=\ln\mathbb{E}[e^{\boldsymbol{\xi}^{T}\boldsymbol{X}}]=\boldsymbol{\xi}^{T}\boldsymbol{\eta}+\frac{1}{2}\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}+\ln[2\Phi(\sqrt{\frac{\pi}{2}}\boldsymbol{\delta}^{T}\operatorname{diag}(\boldsymbol{\sigma})\boldsymbol{\xi}]\,,$
(11)

where
$\boldsymbol{\delta}=\displaystyle\frac{1}{\sqrt{\frac{\pi}{2}(1+\boldsymbol{\alpha}^{T}\boldsymbol{C}\boldsymbol{\alpha})}}\boldsymbol{C}\boldsymbol{\alpha}$,
and
$\boldsymbol{C}=\operatorname{diag}(\boldsymbol{\sigma})^{-1}\boldsymbol{\Sigma}\operatorname{diag}(\boldsymbol{\sigma})^{-1}$.

Since in this framework we work with a centered variable around its mean
$\mkern 1.5mu\overline{\mkern-1.5mu\boldsymbol{X}\mkern-1.5mu}\mkern
1.5mu=\boldsymbol{X}-\boldsymbol{\mu_{X}}$, its CGF becomes

$G_{\mkern 1.5mu\overline{\mkern-1.5mu\boldsymbol{X}\mkern-1.5mu}\mkern
1.5mu}(\boldsymbol{\xi})=-\boldsymbol{\xi}^{T}\boldsymbol{\mu_{X}}+\boldsymbol{\xi}^{T}\boldsymbol{\eta}+\frac{1}{2}\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}+\ln[2\Phi(\sqrt{\frac{\pi}{2}}\boldsymbol{\delta}^{T}\operatorname{diag}(\boldsymbol{\sigma})\boldsymbol{\xi}]\,.$
(12)

Exploiting Expression (9) and denoting
$\boldsymbol{\hat{\mu}}=\boldsymbol{\mu_{X}}-\boldsymbol{\eta}$ and
$\boldsymbol{\xi}=r\boldsymbol{\theta}$, we can write

$G_{\mkern 1.5mu\overline{\mkern-1.5mu\boldsymbol{X}\mkern-1.5mu}\mkern
1.5mu}(r,\boldsymbol{\theta})=-r\boldsymbol{\theta}^{T}\boldsymbol{\hat{\mu}}+\frac{r^{2}}{2}\boldsymbol{\theta}^{T}\boldsymbol{\Sigma}\boldsymbol{\theta}+\ln[2\Phi(\sqrt{\frac{\pi}{2}}r\boldsymbol{\hat{\mu}}^{T}\boldsymbol{\theta})],$
(13)

For small $r$, using the Taylor expansion of
$\ln(1+z)=z-\frac{z^{2}}{2}+\ldots$ and
$2\Phi(z)=1+\mathrm{erf}(\frac{z}{\sqrt{2}})=1+\sqrt{\frac{2}{\pi}}z+\ldots$,
we have
$\ln[2\Phi(\sqrt{\frac{\pi}{2}}r\boldsymbol{\hat{\mu}}^{T}\boldsymbol{\theta})]=\ln(1+\mathrm{erf}(\sqrt{\frac{\pi}{2}}r\boldsymbol{\hat{\mu}}^{T}\boldsymbol{\theta}))\simeq\ln(1+\sqrt{\frac{2}{\pi}}\sqrt{\frac{\pi}{2}}r\boldsymbol{\hat{\mu}}^{T}\boldsymbol{\theta})\simeq
r\boldsymbol{\hat{\mu}}^{T}\boldsymbol{\theta}-\frac{1}{2}r^{2}\boldsymbol{\theta}^{T}\boldsymbol{\hat{\mu}}\boldsymbol{\hat{\mu}}^{T}\boldsymbol{\theta}$.
Therefore, Expression (13) can be approximated as follows

$G_{\mkern 1.5mu\overline{\mkern-1.5mu\boldsymbol{X}\mkern-1.5mu}\mkern
1.5mu}(r,\boldsymbol{\theta})\simeq\frac{r^{2}}{2}\boldsymbol{\theta}^{T}\left(\boldsymbol{\Sigma}-\boldsymbol{\hat{\mu}}\boldsymbol{\hat{\mu}}^{T}\right)\boldsymbol{\theta},$
(14)

where $\boldsymbol{\Sigma}-\boldsymbol{\hat{\mu}}\boldsymbol{\hat{\mu}}^{T}$
is exactly the covariance matrix of the skew-normal distribution
$\boldsymbol{\Sigma_{X}}$ as in (10). Hence, for small $r$, the optimal
solution obtained from Problem (5) collapses to the first principal component,
namely $\boldsymbol{\hat{\theta}_{1}}=PC1$ and
$\boldsymbol{\hat{\theta}_{2}}=-PC1$. Fig. 2 exhibits a numerical example of
the CGF maximization procedure in the case of a bivariate skew-normal random
variable, when $r$ is small.

Figure 2: Bivariate skew-normal random variable: the dashed lines represent
its isoprobability contours, the red solid lines are the isocurves of
$G_{\boldsymbol{X}}(r,\boldsymbol{\theta})=-r\boldsymbol{\mu}^{T}\boldsymbol{\theta}+\frac{r}{2}\boldsymbol{\theta}^{T}\boldsymbol{\Sigma}\boldsymbol{\theta}+\ln[2\Phi(\sqrt{\frac{\pi}{2}}r\boldsymbol{\mu}^{T}\boldsymbol{\theta})]$,
for small $r$.

For large $r$, we can examine the behavior of Eq. (13) using the asymptotic
expansion of the error function (see, e.g., Copson and Copson, 2004) and
distinguishing the case where
$\boldsymbol{\hat{\mu}}^{T}\boldsymbol{\theta}\geq 0$ and
$\boldsymbol{\hat{\mu}}^{T}\boldsymbol{\theta}<0$. For this aim, we consider
the following asymptotic expansion of the complementary error function for
large real $z$

$\mathrm{erfc}(z)=\frac{e^{-z^{2}}}{z{\sqrt{\pi}}}\sum_{n=0}^{\infty}(-1)^{n}{\frac{(2n-1)!!}{\left(2z^{2}\right)^{n}}}=\frac{e^{-z^{2}}}{z{\sqrt{\pi}}}\left(1-\frac{1}{2z^{2}}+\frac{3}{4z^{4}}+\ldots\right)$
(15)

where $(2n-1)!!=(2n-1)\cdot(2n-3)\cdots 3\cdot 1$. Therefore, for
$z\rightarrow+\infty$

$\mathrm{erf}(z)=1-\mathrm{erfc}(z)\simeq
1-\frac{e^{-z^{2}}}{z{\sqrt{\pi}}}+\frac{e^{-z^{2}}}{2z^{3}{\sqrt{\pi}}}$ (16)

This means that

$2\Phi(z)=1+\mathrm{erf}\left(\frac{z}{\sqrt{2}}\right)\xrightarrow[z\to+\infty]{}2$
(17)

Similarly, we obtain

$2\Phi(-z)=1+\mathrm{erf}\left(\frac{-z}{\sqrt{2}}\right)=\mathrm{erfc}\left(\frac{-z}{\sqrt{2}}\right)\xrightarrow[z\to+\infty]{}-\frac{\sqrt{2}e^{-\frac{z^{2}}{2}}}{z{\sqrt{\pi}}}$
(18)
$\ln[2\Phi(-z)]=\ln\left[1+\mathrm{erf}\left(\frac{-z}{\sqrt{2}}\right)\right]=\ln\left[\mathrm{erfc}\left(\frac{-z}{\sqrt{2}}\right)\right]\simeq\ln\left[-\frac{\sqrt{2}e^{-\frac{z^{2}}{2}}}{z{\sqrt{\pi}}}\right]\simeq-\frac{z^{2}}{2}$
(19)

Then, examining Eq. (13) for
$\boldsymbol{\hat{\mu}}^{T}\boldsymbol{\theta}\geq 0$ and large $r$, we have

$G_{\mkern 1.5mu\overline{\mkern-1.5mu\boldsymbol{X}\mkern-1.5mu}\mkern
1.5mu}(r,\boldsymbol{\theta})\simeq\frac{r^{2}}{2}\boldsymbol{\theta}^{T}\boldsymbol{\Sigma}\boldsymbol{\theta}\,,$
(20)

while, for $\boldsymbol{\hat{\mu}}^{T}\boldsymbol{\theta}<0$ and large $r$, we
can write

$G_{\mkern 1.5mu\overline{\mkern-1.5mu\boldsymbol{X}\mkern-1.5mu}\mkern
1.5mu}(r,\boldsymbol{\theta})\simeq-r\boldsymbol{\theta}^{T}\boldsymbol{\hat{\mu}}+\frac{r^{2}}{2}\boldsymbol{\theta}^{T}\boldsymbol{\Sigma}\boldsymbol{\theta}-\frac{1}{2}\frac{\pi}{2}r^{2}\boldsymbol{\theta}^{T}\boldsymbol{\hat{\mu}}\boldsymbol{\hat{\mu}}^{T}\boldsymbol{\theta}\simeq\frac{r^{2}}{2}\boldsymbol{\theta}^{T}(\boldsymbol{\Sigma}-\frac{\pi}{2}\boldsymbol{\hat{\mu}}^{T}\boldsymbol{\hat{\mu}})\boldsymbol{\theta}$
(21)

In Fig. 3, we show a graphical representation of the CGF maximization
procedure in the case of a bivariate skew-normal random variable. The red
solid lines are the isocurves of $G_{\boldsymbol{X}}(r,\boldsymbol{\theta})$,
and the dashed lines represent the isoprobability contours of
$\boldsymbol{X}\sim SN(\boldsymbol{0},\boldsymbol{\Sigma})$, where the
diagonal entries of $\boldsymbol{\Sigma}$ are 1.2 and 0.5143, while the other
entries are equal to 0, and $\boldsymbol{\alpha}=(4.365,-1.455)$.

Figure 3: Bivariate skew-normal random variable: the dashed lines represent
its isoprobability contours, the red solid lines are the isocurves of
$G_{\boldsymbol{X}}(r,\boldsymbol{\theta})=-r\boldsymbol{\mu}^{T}\boldsymbol{\theta}+\frac{r}{2}\boldsymbol{\theta}^{T}\boldsymbol{\Sigma}\boldsymbol{\theta}+\ln[2\Phi(\sqrt{\frac{\pi}{2}}r\boldsymbol{\mu}^{T}\boldsymbol{\theta})]$,
for large $r$.

###  2.3 Maximizing the non-parametric Cumulant Generating Function

In this section we address Problem (5) in the case of a generic discrete
multivariate random variable $\boldsymbol{X}_{t}=(X_{1,t},\ldots,X_{n,t})$
defined on a discrete state space, where we assume $T$ states of nature, each
with probability $\pi_{t}$ with $t=1,\ldots,T$. Therefore, to find the optimal
direction $\widehat{\boldsymbol{\theta}}$ that maximizes the non-parametric
Cumulant Generating Function, we solve the following optimization problem

$\begin{array}[]{lll}\displaystyle\max_{\boldsymbol{\theta}}&\displaystyle
G_{\boldsymbol{X}}(r,\boldsymbol{\theta})=\ln\biggl{(}\sum_{t=1}^{T}\pi_{t}e^{r\boldsymbol{\theta}^{T}\boldsymbol{X}_{t}}\biggr{)}\\\
\mbox{s.t.}&&\\\ &\boldsymbol{\theta}^{T}\boldsymbol{\theta}=1&\end{array}$
(22)

First, in the following theorem we show that the objective function of this
problem is convex. Let us denote, for the sake of notation,

$g(\boldsymbol{\xi})=G_{\boldsymbol{X}}(\boldsymbol{\xi})=\ln\biggl{(}\sum_{t=1}^{T}\pi_{t}e^{\boldsymbol{\xi}^{T}\boldsymbol{X}_{t}}\biggr{)}=G_{\boldsymbol{X}}(r,\boldsymbol{\theta})=\ln\biggl{(}\sum_{t=1}^{T}\pi_{t}e^{r\boldsymbol{\theta}^{T}\boldsymbol{X}_{t}}\biggr{)}\,,$
(23)

where $\boldsymbol{\xi}=r\boldsymbol{\theta}$.

######  Theorem 1.

Let $g(\boldsymbol{\xi})$ be as in (23), where
$\boldsymbol{\xi}\in\mathbb{R}^{n}$. Then, $g(\boldsymbol{\xi})$ is a convex
function.

###### Proof.

Let $\boldsymbol{\xi},\boldsymbol{\zeta}\in\mathbb{R}^{n}$, and
$\gamma\in[0,1]$. Thus,

$g(\gamma\boldsymbol{\xi}+(1-\gamma)\boldsymbol{\zeta})=\ln\biggl{(}\sum_{t=1}^{T}\pi_{t}e^{(\gamma\boldsymbol{\xi}^{T}+(1-\gamma)\boldsymbol{\zeta}^{T})\boldsymbol{X}_{t}}\biggr{)}=\ln\biggl{(}\sum_{t=1}^{T}\pi_{t}e^{\gamma\boldsymbol{\xi}^{T}\boldsymbol{X}_{t}}\,\,e^{(1-\gamma)\boldsymbol{\zeta}^{T}\boldsymbol{X}_{t}}\biggr{)}$

Now, let $u_{t}=\pi_{t}e^{\boldsymbol{\xi}^{T}\boldsymbol{X}_{t}}$ and
$v_{t}=\pi_{t}e^{\boldsymbol{\zeta}^{T}\boldsymbol{X}_{t}}$ $\forall
t=1,\ldots,T$, hence

$g(\gamma\boldsymbol{\xi}+(1-\gamma)\boldsymbol{\zeta})=\ln\biggl{(}\sum_{t=1}^{T}u_{t}^{\gamma}v_{t}^{1-\gamma}\biggr{)}$
(24)

From Hölder’s inequality, we have

$\sum_{t=1}^{T}\lvert x_{t}y_{t}\rvert\leq\biggl{(}\sum_{t=1}^{T}\lvert
x_{t}\rvert^{p}\biggr{)}^{\frac{1}{p}}\biggl{(}\sum_{t=1}^{T}\lvert
y_{t}\rvert^{q}\biggr{)}^{\frac{1}{q}}\,,$

where $\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^{T}$ and $p,q\in[1,\infty)$,
with $\displaystyle\frac{1}{p}+\frac{1}{q}=1$.

Now, considering $x_{t}=u_{t}^{\gamma}$, $y_{t}=v_{t}^{1-\gamma}$,
$\displaystyle\frac{1}{p}=\gamma$ and $\displaystyle\frac{1}{q}=1-\gamma$, we
can apply Hölder’s inequality to Expression (24) as follows

$\begin{split}g(\gamma\boldsymbol{\xi}+(1-\gamma)\boldsymbol{\zeta})=\ln\biggl{(}\sum_{t=1}^{T}u_{t}^{\gamma}v_{t}^{1-\gamma}\biggr{)}&\leq\ln\biggl{(}\biggl{[}\sum_{t=1}^{T}u_{t}^{\gamma\frac{1}{\gamma}}\biggr{]}^{\gamma}\biggl{[}\sum_{t=1}^{T}v_{t}^{1-\gamma\frac{1}{1-\gamma}}\biggr{]}^{1-\gamma}\biggr{)}=\\\
&=\ln\biggl{(}\biggl{[}\sum_{t=1}^{T}u_{t}\biggr{]}^{\gamma}\biggl{[}\sum_{t=1}^{T}v_{t}\biggr{]}^{1-\gamma}\biggr{)}=\\\
&=\gamma\ln\biggl{(}\sum_{t=1}^{T}u_{t}\biggr{)}+(1-\gamma)\ln\biggl{(}\sum_{t=1}^{T}v_{t}\biggr{)}=\\\
&=\gamma g(\boldsymbol{\xi})+(1-\gamma)\,g(\boldsymbol{\zeta})\end{split}$

which completes the proof. ∎

Therefore, Problem (22) consists in maximizing a convex function on an
$n$-sphere, that is a nonconvex set. However, thanks to the following
proposition we can relax such an $n$-sphere into an $n$-ball.

######  Proposition 2.

Problem (22) is equivalent to

$\begin{array}[]{lll}\displaystyle\max_{\boldsymbol{\theta}}&\displaystyle
G_{\boldsymbol{X}}(r,\boldsymbol{\theta})=\ln\biggl{(}\sum_{t=1}^{T}\pi_{t}e^{r\boldsymbol{\theta}^{T}\boldsymbol{X}_{t}}\biggr{)}\\\
\mbox{s.t.}&&\\\ &\boldsymbol{\theta}^{T}\boldsymbol{\theta}\leq
1&\end{array}$ (25)

###### Proof.

Since the $n$-ball $\boldsymbol{\theta}^{T}\boldsymbol{\theta}\leq 1$ can be
seen as the convex hull of the $n$-sphere
$\boldsymbol{\theta}^{T}\boldsymbol{\theta}=1$, we can apply Theorem 32.2 of
Rockafellar (1970), which concludes the proof. ∎

Thus, Problem (25) is a concave programming problem, which is NP hard (Benson,
1995). Indeed, the problem of globally maximizing a convex function on a
convex set may have many local minima, hence finding the global maximum is a
computationally difficult problem, and several approaches have been proposed
in the literature to address this problem (see, e.g., Pardalos and Rosen,
1986).

###  2.4 A Heuristic for maximizing CGF

Similar to Bernacchia et al. (2008), to solve Problem (25) we use a multistart
method, which consists in generating several random initial points belonging
to the $n$-sphere to determine a point in a neighborhood of a global maximum,
and in using a local maximizer to efficiently determine a global maximum. We
point out that we randomly generate initial points on the $n$-sphere, because,
if there exists at least a global maximum of CGF on an $n$-ball, this belongs
to its frontier (see Corollary 32.3.1 of Rockafellar, 1970). Furthermore,
since the constraint of Problem (25) defines a convex set, we use the
projected gradient algorithm as the local maximizer (see Trendafilov and
Jolliffe, 2006). In Table 1 we summarize the CGF maximization procedure.

1. Generate $N$ starting points (directions) $\boldsymbol{\theta}_{0}^{(j)}$ with $j=1,\ldots,N$ belonging to the unit $n$-sphere (see Knuth, 2014). Fix the step size $\delta_{i}=\delta=\frac{1}{r}$ and the tolerance $\epsilon$ sufficiently small (see Nocedal and Wright, 1999). 2. for $j=1,\ldots,N$ 3. Set $i=0$ 4. while $\left\|\boldsymbol{\theta}_{i+1}^{(j)}-\boldsymbol{\theta}_{i}^{(j)}\right\|>\epsilon$ 5. Compute the ascent direction as $\nabla_{\boldsymbol{\theta}_{i}}G_{\boldsymbol{X}}(r,\boldsymbol{\theta}_{i}^{(j)})$ 6. Update $\boldsymbol{\theta}_{i+1}^{(j)}=\boldsymbol{\theta}_{i}^{(j)}+\delta_{i}\nabla_{\boldsymbol{\theta}_{i}}G_{\boldsymbol{X}}(r,\boldsymbol{\theta}_{i}^{(j)})$ 7. Project by rescaling $\boldsymbol{\theta}_{i+1}^{(j)}$ to unit norm, i.e., $\boldsymbol{\theta}_{i+1}^{(j)}=\frac{\boldsymbol{\theta}_{i+1}^{(j)}}{\left\|\boldsymbol{\theta}_{i+1}^{(j)}\right\|}$ 8. Update $i=i+1$ 9. end while 10. end for

Table 1: Pseudocode of the CGF maximization procedure

More precisely, in Expression (23), we assume $\pi_{t}=\frac{1}{T}$ for
$t=1,\cdots,T$, and, therefore, the CGF of a discrete multivariate variable
$\boldsymbol{X}_{t}=(X_{1,t},\ldots,X_{n,t})$ becomes

$G_{\boldsymbol{X}}(r,\boldsymbol{\theta})=\ln\displaystyle\frac{1}{T}\sum_{t=1}^{T}e^{r\boldsymbol{\theta}^{T}\boldsymbol{X}_{t}}\,.$
(26)

For a fixed $r$, we set $N=1000$ and $\delta_{i}=\frac{1}{r}$, where $N$ is
the number of starting points of the multistart heuristic and $\delta_{i}$ is
the step size of the projected gradient algorthm (see Table 1). Furthermore,
from Expression (26), for a given $r$, the gradient of
$G_{\boldsymbol{X}}(r,\boldsymbol{\theta})$ is

$\nabla_{\boldsymbol{\theta}}G_{\boldsymbol{X}}(r,\boldsymbol{\theta})=r\displaystyle\frac{\frac{1}{T}\sum_{t=1}^{T}\boldsymbol{X}_{t}e^{r\boldsymbol{\theta}^{T}\boldsymbol{X}_{t}}}{\frac{1}{T}\sum_{t=1}^{T}e^{r\boldsymbol{\theta}^{T}\boldsymbol{X}_{t}}}$
(27)

Therefore, in Step 6 of the pseudocode in Table 1, the iterative scheme to
find a local optimum is defined by

$\boldsymbol{\theta}_{i+1}-\boldsymbol{\theta}_{i}=\delta_{i}\nabla_{\boldsymbol{\theta}_{i}}G_{\boldsymbol{X}}(r,\boldsymbol{\theta}_{i})=\displaystyle\frac{\sum_{t=1}^{T}\boldsymbol{X}_{t}e^{r\boldsymbol{\theta}_{i}^{T}\boldsymbol{X}_{t}}}{\sum_{t=1}^{T}e^{r\boldsymbol{\theta}_{i}^{T}\boldsymbol{X}_{t}}}\,,$
(28)

while, in Step 7, $\boldsymbol{\theta}_{i+1}$ is normalized to one.

Thus applying the CGF maximization procedure to Problem (25), we can obtain
the local maxima $\boldsymbol{\hat{\theta}}=\boldsymbol{\hat{\theta}}(r)$,
which, except for symmetric distributions (see Section 2.2.1), depend on $r$.
As discussed in Section 2, if $r$ (the distance between data points and the
center) is small, then the CGF maximization procedure essentially picks the
first principal component of the classical PCA. Whereas, if $r$ is large, then
the optimal directions $\boldsymbol{\hat{\theta}}$ maximizing
$G_{\boldsymbol{X}}(r,\boldsymbol{\theta})$ strongly depend on the higher-
order cumulants. On the other hand, higher values of $r$ produce a less
accurate estimate of the function $G(r,\boldsymbol{\theta})$, since the data
sample is finite and a few outliers could heavily influence the higher-order
cumulants. Therefore, when setting the value of $r$, one must consider the
trade-off between finding $\boldsymbol{\hat{\theta}}(r)$ with large $r$ (thus
involving higher-order cumulants in the procedure) and obtaining an accurate
estimate of (26). As suggested by Bernacchia et al. (2008), we set a value
$\bar{r}$ such that the CGF estimation error is limited. We identify such an
error by means of the relative variance of the CGF estimate,
$\varepsilon_{G}^{2}$, that, as shown in Appendix A, assuming i.i.d. normally
distributed random vectors $\boldsymbol{X}_{t}\sim
N_{n}(\boldsymbol{0},\boldsymbol{\Sigma})$ $\forall t=1,\ldots,T$, is as
follows

$\varepsilon_{G}^{2}=\frac{\mathbb{V}\mathrm{ar}[G]}{(\mathbb{E}[G])^{2}}\simeq\frac{4}{T}\frac{e^{\bar{r}^{2}\lambda_{1}}-1}{\bar{r}^{4}(\lambda_{1})^{2}}\,,$
(29)

where $T$ is the number of data points and $\lambda_{1}$ is the largest
eigenvalue computed by standard PCA. In the experimental analysis, $\bar{r}$
is found through (29) by setting $\varepsilon_{G}=10\%$.

##  3 Outlier detection methodologies

In this section, we introduce the proposed outlier detection algorithm for
multivariate data, named the MaxCGF algorithm, consisting in finding outliers
in univariate projections of such data. Our approach is compared with two
other projection-based methods, i.e., one developed by Peña and Prieto (2001,
2007) and the other developed by Domino (2020). Their main difference relies
on the selected directions onto which the data are projected. More precisely,
Peña and Prieto (2001, 2007) consider the directions for which the kurtosis of
the projected data shows the highest and lowest values. Domino (2020) takes
into account the projection directions for which a multivariate series
exhibits the highest absolute value of the 4th order cumulant. We propose a
more general approach that is based on the projection directions that maximize
the Cumulant Generating Function of the multivariate data. In the following we
report the main steps of our procedure.

  1. 1.

Preprocess the original data $\boldsymbol{X}$, yielding centered data
$\boldsymbol{Y}$. Note that while both Peña and Prieto (2001, 2007) and Domino
(2020) standardize $\boldsymbol{X}$, we only center them.

  2. 2.

Find the directions that maximize CGF of $\boldsymbol{Y}$.

  3. 3.

Compute the univariate projection $\boldsymbol{Z}$ of $\boldsymbol{Y}$ in the
directions identified in Step 2.

  4. 4.

Remove outliers. Note that to identify whether a generic element $z_{t}$ of
$\boldsymbol{Z}=\\{z_{t}\\}_{t=1,\ldots,T}$ is an outlier, similar to Peña and
Prieto (2001, 2007) and Domino (2020), we compute the following quantity

$q_{t}=\frac{|z_{t}-\mbox{median}(\boldsymbol{Z})|}{MAD(\boldsymbol{Z})},\quad
t=1,\ldots,T\,,$ (30)

where $T$ is the length of the time series, and $MAD(\boldsymbol{Z})$
represents the median absolute deviation of $\boldsymbol{Z}$. Then, a generic
outcome $z_{t}$ is classified as an outlier if its corresponding $q_{t}$
exceeds a fixed threshold $\beta$. Such a threshold is selected to cover the
whole detection range both in terms of True Positive Rate and in terms of
False Positive Rate (for more details, see Section 4). In order to determine
the optimal $\beta$, we consider equally spaced values of $\beta$ ranging from
0.5 to 10 with a step of 0.25.

  5. 5.

Repeat Steps 2, 3 and 4 until the kurtosis of the univariate projection
$\boldsymbol{Z}$ increases (see Domino, 2020).

In Table 2 we summarize the MaxCGF algorithm (pseudocode) for outlier
detection.

1. Fix the threshold $\beta$; 2. Center data, i.e., $\boldsymbol{Y}^{(0)}=\boldsymbol{X}-\boldsymbol{\mu}$; 3. Find the $N$ directions $\hat{\boldsymbol{\theta}}_{1}^{(0)},\hat{\boldsymbol{\theta}}_{2}^{(0)},\ldots,\hat{\boldsymbol{\theta}}_{N}^{(0)}$, maximizing CGF of $\boldsymbol{Y}^{(0)}$ (see Table 1); 4. for $j=1:N$ (i.e., for each direction $\hat{\boldsymbol{\theta}}_{j}^{(0)}$) 5. Set $i=0$; 6. Project $\boldsymbol{Y}^{(0)}$ on $\hat{\boldsymbol{\theta}}_{j}^{(0)}$, thus obtaining the vector $\boldsymbol{Z}_{j}^{(0)}=\boldsymbol{Y}^{(0)}\hat{\boldsymbol{\theta}}_{j}^{(0)}$; 7. Compute $\mathrm{Kur}_{0}=\mathbb{K}\mathrm{ur}(\boldsymbol{Z}_{j}^{(0)})$, i.e., the kurtosis of $\boldsymbol{Z}_{j}^{(0)}$; 8. while $\mathrm{Kur}_{i}<\mathrm{Kur}_{i-1}$ or $i=0$ 9. Compute the vector $\boldsymbol{q}=\\{q_{t}\\}_{t=1,\ldots,T}$ as in (30) for $\boldsymbol{Z}_{j}^{(i)}$; 10. Remove the outliers when $q_{t}>\beta$, thus obtaining $\boldsymbol{Y}^{(i+1)}$; 11. From $\boldsymbol{Y}^{(i+1)}$, compute $\hat{\boldsymbol{\theta}}_{j}^{(i+1)}$ by using the algorithm in Table 1; 12. Project $\boldsymbol{Y}^{(i+1)}$ on $\hat{\boldsymbol{\theta}}_{j}^{(i+1)}$, thus obtaining the vector $\boldsymbol{Z}_{j}^{(i+1)}=\boldsymbol{Y}^{(i+1)}\hat{\boldsymbol{\theta}}_{j}^{(i+1)}$; 13. Compute $\mathrm{Kur}_{i+1}=\mathbb{K}\mathrm{ur}(\boldsymbol{Z}_{j}^{(i+1)})$, i.e., the kurtosis of $\boldsymbol{Z}_{j}^{(i+1)}$; 14. Update $i=i+1$ 15. end while 16. $\boldsymbol{Y}^{(0)}=\boldsymbol{Y}^{(i+1)}$, $\hat{\boldsymbol{\theta}}_{j}^{(0)}=\hat{\boldsymbol{\theta}}_{j}^{(i+1)}$ 17. end for

Table 2: Pseudocode of the MaxCGF algorithm

In the next section, we compare the outlier detection ability of our method
with that of two alternative methods proposed by Peña and Prieto (2001, 2007)
and Domino (2020).

##  4 Experimental analysis

In this section, we provide a thorough empirical analysis, where the three
outlier detection procedures discussed in Section 3 are tested and compared
both using simulated data (see Section 4.1) and financial real-world data (see
Section 4.2). Performances are evaluated using the Receiver Operating
Characteristic (ROC) curves, which are built by plotting the True Positive
Rate (TPR) vs. the False Positive Rate (FPR) of each methodology for different
levels of the threshold $\beta$. TPR is the rate of truly detected outliers
w.r.t. all possible outliers, and FPR (i.e., the type I error) is the rate of
falsely detected outliers w.r.t. all the non-outlier data. Furthermore, we
also provide two other performance measures related to the ROC curve, namely
the Area Under the Curve (AUC) and Youden’s J Statistic (YJS). AUC is the area
underneath the ROC curve: if a method perfectly distinguishes between outlier
and non-outlier data, then its AUC is equal to 1; if, on the other hand, AUC
is less than or equal to 0.5, then this measure is uninformative, i.e., there
is no difference in performance between the analyzed method and one that
relies on random choices. YJS represents the difference between TPR and FPR,
and again, the higher, the better. The Best Cutoff Value (BCV) is the largest
YJS, and $\beta^{*}$ is the threshold value corresponding to BCV. Although
$\beta^{*}$ is not a performance measure in the strict sense, it provides
information about the practical use of algorithms and allows one to identify
whether there is an optimal range of $\beta$ for which algorithms work best.

Note that in our experiments we report the computational results for all the
three methods only in three instances (standard normal, normal, and real-world
data from the Dow Jones market). In the remaining cases, we report only the
computational results of two methods since the Peña-Prieto algorithm does not
provide reasonable results, and the generated ROC curves do not show a
nondecreasing trend at some points, or yield AUC values less than 0.5.

We implemented all the experiments on a workstation with Intel(R) Xeon(R) CPU
(E5-2623 v4, 2.6 GHz, 64 Gb RAM) under MS Windows 10, using MATLAB 9.12.0.

###  4.1 Computational results for simulated data

Artificial simulated returns are drawn from four different distributions:
standard normal (Section 4.1.1), normal (Section 4.1.2), skew-normal (Section
4.1.3) and Student’s t (Section 4.1.4), with $n=30$ (the number of marginals,
i.e., assets) and $T=500$ (the number of scenarios), which corresponds to 2
financial years. Following Domino (2020), the outliers matrix
$\boldsymbol{O}_{(0.1T\times 0.5n)}$ has dimensions $0.1\ T\times 0.5\ n$.
Below, we report the main steps for generating outliers and substituting them
into the ordinary data.

  1. 1.

Compute the sample covariance matrix $\boldsymbol{\Sigma}$ from a real-world
dataset (here, the Dow Jones dataset from 07/2004 to 07/2006, see Section
4.2).

  2. 2.

Generate the ordinary data $\boldsymbol{X}_{(T\times n)}$ with $T$ scenarios
from an $n$-dimensional random vector using $\boldsymbol{\Sigma}$ (in the
standard normal case, $\boldsymbol{\Sigma}=\boldsymbol{I}$); generate the
outlier data $\boldsymbol{O}_{(0.1T\times 0.5n)}$ using $\boldsymbol{\Sigma}$
multiplied by 15.

  3. 3.

Substitute the elements of $\boldsymbol{O}_{(0.1T\times 0.5n)}$ in
$\boldsymbol{X}_{(T\times n)}$. More precisely, randomly select $n/2$ columns
(with $0.1T$ elements) of $\boldsymbol{X}_{(T\times n)}$ and substitute them
with the columns of $\boldsymbol{O}_{(0.1T\times 0.5n)}$.

####  4.1.1 The standard normal random vector case

For this experimental case, the ordinary and outlier data are
$\boldsymbol{X}_{(T\times n)}\sim N(\boldsymbol{0},\boldsymbol{I})$ and
$\boldsymbol{O}_{(0.1T\times 0.5n)}\sim N(\boldsymbol{0},15\boldsymbol{I})$.

In Table 3, we report the computational results obtained by using the three
methods described in Section 3. More precisely, the Peña-Prieto method obtains
an almost perfect score of 0.9867 both for the AUC and the BCV performance
measures. The Domino algorithm shows the lowest performance w.r.t. the other
two methods, both in terms of AUC ($=0.9316$) and BCV ($=0.7489$). The MaxCGF
method almost achieves the highest values in terms of accuracy and clearly
shows a significant advantage in terms of computational burden. Its running
time is less than 1 minute, compared with 15 minutes and more than 1 hour
spent by the Domino and Peña-Prieto methods, respectively.

Method |  AUC |  BCV |  Time (min.) |  $\boldsymbol{\beta}^{*}$
---|---|---|---|---
MaxCGF |  $0.9843$ |  $0.9533$ |  $0.6032$ |  3.25
Domino |  $0.9316$ |  $0.7489$ |  $21.6774$ |  4.25
Peña-Prieto |  $0.9867$ |  $0.9867$ |  $66.6850$ |  2.25
Table 3: Performances for the standard normal dataset. We mark in green the
best, in yellow the intermediate, and in red the worst results.

Figure 4 reports the ROC curves of the three methods, which, as already noted,
show high outlier detection abilities in this experiment.

Figure 4: ROC curves of the three algorithms for the standard normal market.

####  4.1.2 The normal random vector case

In this case, the ordinary and outlier data are $\boldsymbol{X}_{(T\times
n)}\sim N(\boldsymbol{0},\boldsymbol{\Sigma})$ and
$\boldsymbol{O}_{(0.1T\times 0.5n)}\sim
N(\boldsymbol{0},15\boldsymbol{\Sigma})$.

The empirical results of the ROC analysis are shown in Figure 5 and in Table
4.

Figure 5: ROC curves of the three algorithms for the normal market. Method |  AUC |  BCV |  Time (min.) |  $\boldsymbol{\beta}^{*}$
---|---|---|---|---
MaxCGF |  $0.8811$ |  $0.6067$ |  $0.9592$ |  8.00
Domino |  $0.8809$ |  $0.6311$ |  $14.3293$ |  7.50
Peña-Prieto |  $0.9847$ |  $0.9644$ |  $67.2382$ |  1.75
Table 4: Performances for the normal dataset.

Here, the ranking of the three algorithms is similar to that of the standard
normal dataset. Again, the Peña-Prieto algorithm outperforms the other two
both in terms of AUC and BCV, while the MaxCGF algorithm is the most
efficient; however, the Domino algorithm yields the second best BCV. It is
noteworthy that, compared to the standard normal case, the introduction of a
correlation structure seems to worsen the results for the MaxCGF and Domino
algorithms, while that of Peña-Prieto does not experience noticeable changes.

####  4.1.3 The skew-normal random vector case

For this experiment, the ordinary and outlier data are
$\boldsymbol{X}_{(T\times n)}\sim
SN(\boldsymbol{0},\boldsymbol{\Sigma},\boldsymbol{\alpha})$ and
$\boldsymbol{O}_{(0.1T\times 0.5n)}\sim
SN(\boldsymbol{0},15\boldsymbol{\Sigma},\boldsymbol{\alpha})$, where
$\boldsymbol{\alpha}$ is a vector whose elements are uniformly distributed in
the interval $[-1,4]$.

Method |  AUC |  BCV |  Time (min.) |  $\boldsymbol{\beta}^{*}$
---|---|---|---|---
MaxCGF |  $0.9140$ |  $0.6911$ |  $1.0387$ |  7.25
Domino |  $0.8896$ |  $0.6467$ |  $15.8041$ |  6.50
Peña-Prieto |  - |  - |  - |  -

Table 5: Performances for the skew-normal dataset. Figure 6: ROC curves of the
MaxCGF and Domino algorithms for the skew-normal market.

As shown in Table 5 and in Figure 6, the MaxCGF algorithm provides, in the
case of asymmetric distribution, the highest AUC (=0.9140) and BCV (=0.6467),
and it is the most efficient approach. The Domino method also performs well,
although lower than the MaxCGF method. Its AUC is 0.8896, and its BCV is
0.6467, while its running time is 15 minutes compared with 1 minute spent by
the MaxCGF method.

####  4.1.4 The Student’s t random vector case

Here, the ordinary and outlier data follow an $n$-variate Student’s t
distribution, $\boldsymbol{X}_{(T\times n)}\sim
St(\boldsymbol{0},\boldsymbol{\Sigma},\nu)$ and $\boldsymbol{O}_{(0.1T\times
0.5n)}\sim St(\boldsymbol{0},15\boldsymbol{\Sigma},\nu)$, where $\nu$ is set
to 5, 10, 100, 1000. Clearly, when $\nu$ is sufficiently high, the Student’s t
random vector approaches the normal one.

Before testing the three methods described in Section 3, we examine the
outlier detection procedure by using the projection directions provided by the
CGF maximization and by the classical PCA.

$\nu$ |  Method |  AUC |  BCV |  Time (min.)
---|---|---|---|---
5 |  Classical PCA |  $0.7757$ |  $0.4467$ |  $0.8389$
|  MaxCGF |  $0.8544$ |  $0.5933$ |  $22.5747$
10 |  Classical PCA |  $0.7988$ |  $0.4089$ |  $0.7955$
|  MaxCGF |  $0.8333$ |  $0.5356$ |  $10.0730$
100 |  Classical PCA |  $0.8434$ |  $0.5600$ |  $0.7952$
|  MaxCGF |  $0.8525$ |  $0.6311$ |  $21.8163$
1000 |  Classical PCA |  $0.8504$ |  $0.5756$ |  $0.7755$
|  MaxCGF |  $0.8409$ |  $0.5489$ |  $17.5729$

Table 6: Outlier detection performances using the CGF maximization and the
classical PCA for the Student’s t random vector.

As shown in Table 6, the latter is much more efficient in terms of running
time, but, except for $\nu=1000$, the MaxCGF method is more precise, as shown
by the values of AUC and BCV. Furthermore, the accuracy of the classical PCA
tend to increase with $\nu$, namely when the simulated data tend to have a
normal distribution. Figure 7 exhibits the ROC curves of such experiments.

(a) $\nu=5$

(b) $\nu=10$

(c) $\nu=100$

(d) $\nu=1000$

Figure 7: ROC curves using the CGF maximization and the classical PCA for the
Student’s t random vector.

For these experiments, we also compare the performances of the Peña-Prieto,
Domino, and MaxCGF methods considering the degrees of freedom of the Student’s
t random vector, $\nu=10$ and $\nu=30$. Note that, as mentioned at the
beginning of Section 4, the Peña-Prieto algorithm does not provide reasonable
results.

Method |  AUC |  BCV |  Time (min.) |  $\boldsymbol{\beta}^{*}$
---|---|---|---|---
MaxCGF |  $0.8333$ |  $0.5356$ |  $10.0730$ |  7.25
Domino |  $0.8498$ |  $0.5622$ |  $8.6374$ |  7.50
Peña-Prieto |  - |  - |  - |  -

Table 7: Performances for the Student’s t dataset for $\nu=10$.

Method |  AUC |  BCV |  Time (min.) |  $\boldsymbol{\beta}^{*}$
---|---|---|---|---
MaxCGF |  $0.9116$ |  $0.7044$ |  $8.0739$ |  7.75
Domino |  $0.8903$ |  $0.6267$ |  $19.1567$ |  7.25
Peña-Prieto |  - |  - |  - |  -

Table 8: Performances for the Student’s t dataset for $\nu=30$.

Tables 7 and 8 show the values of AUC, BCV, and CPU time for $\nu=10$ and
$\nu=30$, respectively. In the former case, the Domino approach slightly
outperforms the MaxCGF one, while in the latter case, the MaxCGF method yields
better results in terms of AUC, BCV, and running time. Figure 8 exhibits the
two methods’ ROC curves, highlighting elevated outlier detection abilities in
these experiments.

(a) $\nu=10$.

(b) $\nu=30$.

Figure 8: ROC curves of the MaxCGF and Domino algorithms for the Student’s t
random vector.

###  4.2 Computational results for financial real-world data

In this section, we use the three outlier detection methods to identify
financial crises, which, as mentioned in the introduction, can be seen as
periods when markets experience atypical behavior. More precisely, we apply
the outlier detection analysis to six real-world financial datasets belonging
to major stock markets across the world. In Table 9, we provide some details
about these datasets, which consist of daily prices, adjusted for dividends
and stock splits, obtained from _Refinitiv_.

Market | Abbrev. | # assets | Country | From-To
---|---|---|---|---
Dow Jones Industrial Average | DJIA | 30 | USA | 09/04/2019-23/03/2020
Euro Stoxx 50 | STOXX50 | 47 | Eurozone
DAX 30 | DAX | 28 | Germany
CAC 40 | CAC | 40 | France
FTSE 100 | FTSE | 98 | UK
EuroNext 100 | N100 | 98 | Eurozone

Table 9: List of the daily datasets analyzed.

From prices, we use both linear and logarithmic returns for the empirical
analysis. Since the results obtained are practically identical, we report here
only those obtained by means of the linear returns. The data analyzed are
daily returns from April 9, 2019, to March 23, 2020 (approximately one
financial year). We have chosen this time horizon because it contains a period
of high instability due to the COVID-19 pandemic, the effects of which occur
approximately in early February 2020 (see, e.g. Shu et al. (2021)).

Method |  AUC |  BCV |  Time (min.) |  $\boldsymbol{\beta}^{*}$
---|---|---|---|---
MaxCGF |  $0.9057$ |  $0.7156$ |  $5.9174$ |  6.75
Domino |  $0.8282$ |  $0.6178$ |  $25.6869$ |  5.75
Peña-Prieto |  0.7873 |  $0.5511$ |  $84.0113$ |  6.50
Table 10: Performances for DJIA. Method |  AUC |  BCV |  Time (min.) |  $\boldsymbol{\beta}^{*}$
---|---|---|---|---
MaxCGF |  $0.8312$ |  $0.5666$ |  $3.6558$ |  8.00
Domino |  $0.7480$ |  $0.4713$ |  $92.4653$ |  7.50
Peña-Prieto |  - |  - |  - |  -
Table 11: Performances for STOXX50. Method |  AUC |  BCV |  Time (min.) |  $\boldsymbol{\beta}^{*}$
---|---|---|---|---
MaxCGF |  $0.8866$ |  $0.6555$ |  $3.8149$ |  5.75
Domino |  $0.7924$ |  $0.5701$ |  $17.2033$ |  7.00
Peña-Prieto |  - |  - |  - |  -
Table 12: Performances for DAX. Method |  AUC |  BCV |  Time (min.) |  $\boldsymbol{\beta}^{*}$
---|---|---|---|---
MaxCGF |  $0.8326$ |  $0.5187$ |  $3.7390$ |  8.25
Domino |  $0.7391$ |  $0.4637$ |  $75.1684$ |  6.50
Peña-Prieto |  - |  - |  - |  -
Table 13: Performances for CAC. Method |  AUC |  BCV |  Time (min.) |  $\boldsymbol{\beta}^{*}$
---|---|---|---|---
MaxCGF |  $0.8021$ |  $0.5207$ |  $5.4309$ |  7.50
Domino |  $0.7061$ |  $0.4641$ |  $2644.9669$ |  7.50
Peña-Prieto |  - |  - |  - |  -
Table 14: Performances for FTSE. Method |  AUC |  BCV |  Time (min.) |  $\boldsymbol{\beta}^{*}$
---|---|---|---|---
MaxCGF |  $0.8643$ |  $0.5474$ |  $6.2306$ |  7.75
Domino |  $0.6794$ |  $0.3624$ |  $2986.0548$ |  7.75
Peña-Prieto |  - |  - |  - |  -
Table 15: Performances for N100.

In Tables 10-15, we report the computational results for all the datasets
listed in Table 9. These tables show that AUC obtained by the MaxCGF approach
is between 0.8021 and 0.9057, and is always higher than that achieved by the
other two methods.

(a) Value of the Dow Jones index. The red line indicates the beginning of the
crisis.

(b) ROC curves of the three algorithms for the DJIA dataset.

Figure 9: DJIA dataset.

(a) Value of the Eurostoxx50 index. The red line indicates the beginning of
the crisis.

(b) ROC curves of the MaxCGF and Domino algorithms for the STOXX50 dataset.

Figure 10: STOXX50 dataset.

(a) Value of the DAX index. The red line indicates the beginning of the
crisis.

(b) ROC curves of the MaxCGF and Domino algorithms for the DAX dataset.

Figure 11: DAX dataset.

(a) Value of the CAC40 index. The red line indicates the beginning of the
crisis.

(b) ROC curves of the MaxCGF and Domino algorithms for the CAC40 dataset.

Figure 12: CAC40 dataset.

(a) Value of the FTSE100 index. The red line indicates the beginning of the
crisis.

(b) ROC curves of the MaxCGF and Domino algorithms for the FTSE100 dataset.

Figure 13: FTSE100 dataset.

(a) Value of the Euro Next 100 index. The red line indicates the beginning of
the crisis.

(b) ROC curves of the MaxCGF and Domino algorithms for the N100 dataset.

Figure 14: N100 dataset.

In Figures 9-14, we also provide the ROC curves for all the financial market
analyzed. Interestingly, TPR of the MaxCGF algorithm is high for values of FPR
slightly higher than 0.1. This means that our algorithm is able to correctly
detect a large portion of outliers misclassifying only few ordinary data (as
outlier) even for high values of the threshold $\beta$. Furthermore, the
MaxCGF method is always the least time-consuming, only taking a few minutes to
complete the analysis for all the financial datasets analyzed. Conversely,
both the Domino and the Peña-Prieto algorithms are heavily influenced by the
number of assets $n$. Indeed, on the one hand, the Peña-Prieto procedure is
able to produce results only for DJIA. On the other hand, for the Domino
approach, the running time ranges from about 17 minutes for DAX ($n=28$) to
around two days for the N100 ($n=98$). Summing up, as highlighted in the
tables, the MaxCGF method always shows better results for all the examined
performance measures, followed by the Domino method.

##  5 Conclusions

In this paper, we have proposed a non-parametric approach to detect anomalies
in data by examining their univariate projections on appropriate directions
that depend on cumulants of any order. Such directions, where the original
data have been projected, are those that maximize the cumulant generating
function (CGF).

In this respect, we have first refined some theoretical results of Bernacchia
and Naveau (2008) and Bernacchia et al. (2008) investigating the directions
that maximize CGF of data with normal and skew-normal distributions. Then, we
have proved in the general non-parametric case that CGF is a convex function
and characterized the CGF maximization problem on the unit $n$-circle as a
concave minimization problem. Furthermore, we have extended the outlier
detection methodology based on the projections of multivariate data on the
directions obtained by the classical PCA technique to the directions that
maximize CGF. Finally, we have presented an extensive empirical analysis
testing the performance of our outlier detection procedure, named MaxCGF, and
comparing it with two other methods, proposed by Domino (2020) and Peña and
Prieto (2001). From the computational results, we have observed that the Peña-
Prieto method shows high performance for standard normal and normal simulated
data, while, in other cases, its ability to detect outliers seems to be low.
The Domino method typically provides intermediate performance, but favorable
in the case of Student’s t distributed data. The MaxCGF approach also performs
well for normal and Student’s random vectors. However, for skew-normal data or
real-world financial data, the MaxCGF method always shows better results for
all performance measures examined, including efficiency.

##  Appendix A Relative variance of the CGF estimator

The sample estimate of the Cumulant Generating Function (CGF) of a discrete
multivariate variable $\boldsymbol{X}_{t}=(X_{1,t},\ldots,X_{n,t})$ with
$t=1,\ldots,T$ is

$G_{\boldsymbol{X}}(\boldsymbol{\xi})=\ln\displaystyle\frac{1}{T}\sum_{t=1}^{T}e^{\boldsymbol{\xi}^{T}\boldsymbol{X}_{t}}=\ln
M_{\boldsymbol{X}}(\boldsymbol{\xi})\,,$

where, therefore, $M_{\boldsymbol{X}}(\boldsymbol{\xi})$ denotes the sample
estimate of the moment generating function. The relative variance of the CGF
estimator is then defined as

$\varepsilon_{G}^{2}=\frac{\mathbb{V}\mathrm{ar}[G]}{(\mathbb{E}[G])^{2}}\,.$
(31)

To find an explicit expression for (31), we exploit the Taylor expansion for
moments of functions of random variables (see, e.g., Benaroya et al., 2005).
More precisely, we use the Taylor expansion around $\mu_{M}=\mathbb{E}[M]$,
namely

$G=\ln
M=\ln\mu_{M}+\frac{1}{\mu_{M}}(M-\mu_{M})-\frac{1}{2}\frac{1}{\mu_{M}^{2}}(M-\mu_{M})^{2}+\ldots$
(32)

Using the first order approximation, we obtain

$\mathbb{E}[G]=\mathbb{E}[\ln
M]\simeq\mathbb{E}[\ln\mu_{M}+\frac{1}{\mu_{M}}(M-\mu_{M})]=\ln\mu_{M}$ (33)
$\mathbb{V}\mathrm{ar}[G]=\mathbb{V}\mathrm{ar}[\ln
M]\simeq\mathbb{V}\mathrm{ar}[\ln\mu_{M}+\frac{1}{\mu_{M}}(M-\mu_{M})]=\frac{1}{\mu_{M}^{2}}\sigma_{M}^{2}\,.$
(34)

Thus,

$\varepsilon_{G}^{2}\simeq\frac{\sigma_{M}^{2}}{\mu_{M}^{2}(\ln\mu_{M})^{2}}\,,$
(35)

where $\mu_{M}=\mathbb{E}[M]$ and $\sigma_{M}^{2}=\mathbb{V}\mathrm{ar}[M]$.

The sample estimate of the moment generating function of a discrete
multivariate variable $\boldsymbol{X}_{t}=(X_{1,t},\ldots,X_{n,t})$ with
$t=1,\ldots,T$ is

$M_{\boldsymbol{X}}(\boldsymbol{\xi})=\displaystyle\frac{1}{T}\sum_{t=1}^{T}e^{\boldsymbol{\xi}^{T}\boldsymbol{X}_{t}}\,,$
(36)

where, as mentioned in Section 2.4, to explicitly find (35), we assume i.i.d.
normally distributed random vectors $\boldsymbol{X}_{t}\sim
N_{n}(\boldsymbol{0},\boldsymbol{\Sigma})$ $\forall t=1,\ldots,T$. Then, the
expectation of the estimator (36) is

$\mathbb{E}[M]=\displaystyle\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}[e^{\boldsymbol{\xi}^{T}\boldsymbol{X}_{t}}]\,,$
(37)

where $\mathbb{E}[e^{\boldsymbol{\xi}^{T}\boldsymbol{X}_{t}}]=\int\cdots\int
e^{\boldsymbol{\xi}\boldsymbol{x}}f_{\boldsymbol{X}_{t}}(\boldsymbol{x})d\boldsymbol{x}$
and
$f_{\boldsymbol{X}_{t}}(\boldsymbol{x})=2\pi^{-\frac{n}{2}}\det(\boldsymbol{\Sigma})^{-\frac{1}{2}}e^{-\frac{1}{2}\boldsymbol{x}^{T}\boldsymbol{\Sigma}^{-1}\boldsymbol{x}}$.
Hence,

$\displaystyle\mathbb{E}[e^{\boldsymbol{\xi}^{T}\boldsymbol{X}_{t}}]$
$\displaystyle=$ $\displaystyle
2\pi^{-\frac{n}{2}}\det(\boldsymbol{\Sigma})^{-\frac{1}{2}}\int\cdots\int
e^{\boldsymbol{\xi}\boldsymbol{x}}e^{-\frac{1}{2}\boldsymbol{x}^{T}\boldsymbol{\Sigma}^{-1}\boldsymbol{x}}d\boldsymbol{x}$
(38) $\displaystyle=$ $\displaystyle
e^{\frac{1}{2}\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}\,2\pi^{-\frac{n}{2}}\det(\boldsymbol{\Sigma})^{-\frac{1}{2}}\int\cdots\int
e^{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\Sigma}\boldsymbol{\xi})^{T}\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\Sigma}\boldsymbol{\xi})}d\boldsymbol{x}$
$\displaystyle=$ $\displaystyle
e^{\frac{1}{2}\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}\,,$

since
$2\pi^{-\frac{n}{2}}\det(\boldsymbol{\Sigma})^{-\frac{1}{2}}\int\cdots\int
e^{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\Sigma}\boldsymbol{\xi})^{T}\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\Sigma}\boldsymbol{\xi})}d\boldsymbol{x}=1$.
Thus, using (38) in (37), we obtain

$\mu_{M}=\mathbb{E}[M]=\displaystyle\frac{1}{T}\,T\,e^{\frac{1}{2}\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}=e^{\frac{1}{2}\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}$
(39)

For the variance of the estimator (36), we have

$\mathbb{V}\mathrm{ar}[M]=\mathbb{E}[M^{2}]-\mathbb{E}[M]^{2}\,,$ (40)

where

$\displaystyle\mathbb{E}[M^{2}]$ $\displaystyle=$
$\displaystyle\mathbb{E}\biggl{[}\frac{1}{T}\sum_{t=1}^{T}e^{\boldsymbol{\xi}^{T}\boldsymbol{X}_{t}}\,\frac{1}{T}\sum_{j=1}^{T}e^{\boldsymbol{\xi}^{T}\boldsymbol{X}_{j}}\biggr{]}$
$\displaystyle=$
$\displaystyle\displaystyle\frac{1}{T^{2}}\mathbb{E}\biggl{[}\sum_{t=1}^{T}e^{2\boldsymbol{\xi}^{T}\boldsymbol{X}_{t}}+\sum_{t\neq
j}^{T}(e^{\boldsymbol{\xi}^{T}\boldsymbol{X}_{t}})(e^{\boldsymbol{\xi}^{T}\boldsymbol{X}_{j}})\biggr{]}$
$\displaystyle=$
$\displaystyle\displaystyle\frac{1}{T^{2}}\sum_{t=1}^{T}\mathbb{E}\biggl{[}e^{2\boldsymbol{\xi}^{T}\boldsymbol{X}_{t}}\biggr{]}+\sum_{t\neq
j}^{T}\mathbb{E}\biggl{[}(e^{\boldsymbol{\xi}^{T}\boldsymbol{X}_{t}})(e^{\boldsymbol{\xi}^{T}\boldsymbol{X}_{j}})\biggr{]}\,.$

Now, similarly to (38), we obtain that

$\mathbb{E}[e^{2\boldsymbol{\xi}\boldsymbol{X}_{t}}]=e^{2\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}\,.$
(41)

Furthermore, since, by assumption, $\boldsymbol{X}_{t}$ are i.i.d. $\forall
t=1,\ldots,T$, we have

$\displaystyle\mathbb{E}[M^{2}]$ $\displaystyle=$
$\displaystyle\displaystyle\frac{1}{T^{2}}(T\mathbb{E}[e^{2\boldsymbol{\xi}^{T}\boldsymbol{X}}]+(T^{2}-T)(\mathbb{E}[e^{\boldsymbol{\xi}^{T}\boldsymbol{X}}])^{2}$
(42) $\displaystyle=$
$\displaystyle\frac{1}{T}\mathbb{E}[e^{2\boldsymbol{\xi}^{T}\boldsymbol{X}}]+(1-\frac{1}{T})(\mathbb{E}[e^{\boldsymbol{\xi}^{T}\boldsymbol{X}}])^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{T}e^{2\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}+(1-\frac{1}{T})e^{\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}\,.$

Substituting (42) in Expression (40), we can write

$\displaystyle\sigma_{M}^{2}=\mathbb{V}\mathrm{ar}[M]$ $\displaystyle=$
$\displaystyle\frac{1}{T}e^{2\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}+(1-\frac{1}{T})e^{\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}-e^{\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}$
(43) $\displaystyle=$
$\displaystyle\frac{1}{T}\left(e^{2\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}-e^{\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}\right)$

Thus, using (39) and (43) in (31),

$\varepsilon_{G}^{2}\simeq\frac{\sigma_{M}^{2}}{\mu_{M}^{2}(\ln\mu_{M})^{2}}=\frac{\frac{1}{T}\left(e^{2\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}-e^{\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}\right)}{e^{\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}(\frac{1}{2}\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi})^{2}}=\frac{1}{T}\frac{e^{\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi}}-1}{(\frac{1}{2}\boldsymbol{\xi}^{T}\boldsymbol{\Sigma}\boldsymbol{\xi})^{2}}\,.$
(44)

Denoting $\boldsymbol{\xi}=r\boldsymbol{\theta}$, we obtain

$\varepsilon_{G}^{2}\simeq\frac{4}{T}\frac{e^{r^{2}\boldsymbol{\theta}^{T}\boldsymbol{\Sigma}\boldsymbol{\theta}}-1}{r^{4}(\boldsymbol{\theta}^{T}\boldsymbol{\Sigma}\boldsymbol{\theta})^{2}}\,.$
(45)

Simplifying the issue, we substitute
$\boldsymbol{\theta}^{T}\boldsymbol{\Sigma}\boldsymbol{\theta}$ with its
largest eigenvalue $\lambda_{1}$ (computed by the standard PCA), and,
therefore, Expression (45) becomes

$\varepsilon_{G}^{2}\simeq\frac{4}{T}\frac{e^{r^{2}\lambda_{1}}-1}{r^{4}(\lambda_{1})^{2}}\,.$
(46)

## References

  * Ané et al. (2008) Ané, T., Ureche-Rangau, L., Gambet, J.-B., and Bouverot, J. (2008).  Robust outlier detection for asia–pacific stock index returns.  Journal of International Financial Markets, Institutions and Money, 18(4):326–343.
  * Arellano-Valle and Azzalini (2008) Arellano-Valle, R. B. and Azzalini, A. (2008).  The centred parametrization for the multivariate skew-normal distribution.  Journal of Multivariate Analysis, 99(7):1362–1382.
  * Azzalini and Capitanio (1999) Azzalini, A. and Capitanio, A. (1999).  Statistical applications of the multivariate skew normal distribution.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61(3):579–602.
  * Azzalini and Valle (1996) Azzalini, A. and Valle, A. D. (1996).  The multivariate skew-normal distribution.  Biometrika, 83(4):715–726.
  * Benaroya et al. (2005) Benaroya, H., Han, S. M., and Nagurka, M. (2005).  Probability models in engineering and science, volume 192.  CRC press.
  * Benson (1995) Benson, H. P. (1995).  Concave minimization: theory, applications and algorithms.  In Handbook of global optimization, pages 43–148. Springer.
  * Bernacchia and Naveau (2008) Bernacchia, A. and Naveau, P. (2008).  Detecting spatial patterns with the cumulant function–part 1: The theory.  Nonlinear Processes in Geophysics, 15(1):159–167.
  * Bernacchia et al. (2008) Bernacchia, A., Naveau, P., Vrac, M., and Yiou, P. (2008).  Detecting spatial patterns with the cumulant function–part 2: An application to el nino.  Nonlinear Processes in Geophysics, 15(1):169–177.
  * Cesarone et al. (2020) Cesarone, F., Mango, F., Mottura, C. D., Ricci, J. M., and Tardella, F. (2020).  On the stability of portfolio selection models.  Journal of Empirical Finance, 59:210–234.
  * Copson and Copson (2004) Copson, E. T. and Copson, E. T. (2004).  Asymptotic expansions.  Number 55. Cambridge university press.
  * Das and Sinha (1986) Das, R. and Sinha, B. K. (1986).  Detection of multivariate outliers with dispersion slippage in elliptically symmetric distributions.  The Annals of Statistics, pages 1619–1624.
  * Domino (2020) Domino, K. (2020).  Multivariate cumulants in outlier detection for financial data analysis.  Physica A: Statistical Mechanics and its Applications, page 124995\.
  * Donoho (1982) Donoho, D. L. (1982).  Breakdown properties of multivariate location estimators.  Technical report, Technical report, Harvard University, Boston. URL http://www-stat. stanford ….
  * Ferguson (1961) Ferguson, T. S. (1961).  On the rejection of outliers.  In Proceedings of the fourth Berkeley symposium on mathematical statistics and probability, volume 1, pages 253–287. University of California Press Berkeley and Los Angeles.
  * Giacometti et al. (2021) Giacometti, R., Torri, G., and Paterlini, S. (2021).  Tail risks in large portfolio selection: penalized quantile and expectile minimum deviation models.  Quantitative Finance, 21(2):243–261.
  * Gnanadesikan and Kettenring (1972) Gnanadesikan, R. and Kettenring, J. R. (1972).  Robust estimates, residuals, and outlier detection with multiresponse data.  Biometrics, pages 81–124.
  * Knuth (2014) Knuth, D. E. (2014).  Art of computer programming, volume 2: Seminumerical algorithms.  Addison-Wesley Professional.
  * Kondor et al. (2007) Kondor, I., Pafka, S., and Nagy, G. (2007).  Noise sensitivity of portfolio selection under various risk measures.  Journal of Banking & Finance, 31(5):1545–1573.
  * Maronna (1976) Maronna, R. A. (1976).  Robust M-estimators of multivariate location and scatter.  The Annals of Statistics, pages 51–67.
  * Meng et al. (2019) Meng, F., Yuan, G., Lv, S., Wang, Z., and Xia, S. (2019).  An overview on trajectory outlier detection.  Artificial Intelligence Review, 52(4):2437–2456.
  * Nocedal and Wright (1999) Nocedal, J. and Wright, S. J. (1999).  Numerical optimization.  Springer.
  * Pardalos and Rosen (1986) Pardalos, P. M. and Rosen, J. B. (1986).  Methods for global concave minimization: A bibliographic survey.  Siam Review, 28(3):367–379.
  * Peña and Prieto (2001) Peña, D. and Prieto, F. J. (2001).  Multivariate outlier detection and robust covariance matrix estimation.  Technometrics, 43(3):286–310.
  * Peña and Prieto (2007) Peña, D. and Prieto, F. J. (2007).  Combining random and specific directions for outlier detection and robust estimation in high-dimensional multivariate data.  Journal of Computational and Graphical Statistics, 16(1):228–254.
  * Reed and Yu (1990) Reed, I. S. and Yu, X. (1990).  Adaptive multiple-band cfar detection of an optical pattern with unknown spectral distribution.  IEEE transactions on acoustics, speech, and signal processing, 38(10):1760–1770.
  * Rockafellar (1970) Rockafellar, R. T. (1970).  Convex Analysis, volume 36.  Princeton University Press.
  * Rousseeuw (1984) Rousseeuw, P. J. (1984).  Least median of squares regression.  Journal of the American Statistical Association, 79(388):871–880.
  * Rousseeuw (1985) Rousseeuw, P. J. (1985).  Multivariate estimation with high breakdown point.  Mathematical Statistics and Applications, 8(283-297):37.
  * Schwager and Margolin (1982) Schwager, S. J. and Margolin, B. H. (1982).  Detection of multivariate normal outliers.  The Annals of Statistics, 10(3):943–954.
  * Shu et al. (2021) Shu, M., Song, R., and Zhu, W. (2021).  The “COVID” crash of the 2020 US Stock market.  The North American Journal of Economics and Finance, 58:101497.
  * Singh and Upadhyaya (2012) Singh, K. and Upadhyaya, S. (2012).  Outlier detection: applications and techniques.  International Journal of Computer Science Issues (IJCSI), 9(1):307.
  * Sinha (1984) Sinha, B. K. (1984).  Detection of multivariate outliers in elliptically symmetric distributions.  The Annals of Statistics, pages 1558–1565.
  * Stahel (1981) Stahel, W. A. (1981).  Breakdown of covariance estimators.  Fachgruppe für Statistik, Eidgenössische Techn. Hochsch.
  * Trendafilov and Jolliffe (2006) Trendafilov, N. T. and Jolliffe, I. T. (2006).  Projected gradient approach to the numerical solution of the scotlass.  Computational Statistics & Data Analysis, 50(1):242–253.
  * Wilks (1963) Wilks, S. S. (1963).  Multivariate statistical outliers.  Sankhyā: The Indian Journal of Statistics, Series A, pages 407–426.