Mean-Variance-Skewness-Kurtosis Analysis

• Mean-variance-skewness-kurtosis problems model asset returns by jointly considering the first four moments to capture asymmetry and heavy tails.
• They employ diverse formulations such as weighted-sum scalarization, feasibility analysis, and cumulant ratios to address nonconvex, quartic optimization challenges.
• Advanced methods like successive convex approximation, DC programming, and tensor-free techniques enable scalable solutions for complex, large-scale moment constraints.

Searching arXiv for the cited MVSK and related papers to ground the article in the current arXiv record. Tool unavailable in this interface, so proceeding with the arXiv records explicitly provided in the data block and citing those papers directly. The mean-variance-skewness-kurtosis problem concerns modeling, inference, or optimization in which the first four moments or cumulants of a random quantity are treated jointly. In portfolio selection, it usually denotes the design of portfolios that reward mean and skewness while penalize variance and kurtosis, motivated by the fact that Gaussian mean-variance descriptions are inadequate when returns are asymmetric and heavy-tailed. In related statistical and physical settings, the same quartet appears as raw central moments, standardized moments, or cumulant ratios, so the problem is simultaneously an optimization problem, a feasibility problem for admissible moment tuples, and an inferential problem for estimating higher-order structure (Zhou et al., 2020, Collaboration et al., 2017, Meer et al., 2023).

1. Definitions and moment conventions

A first source of technical variation is that different literatures use different versions of the four moments. In high-order portfolio optimization, the portfolio return is $\mathbf{w}^\top \mathbf{r}$, with mean, variance, skewness, and kurtosis defined as

$$\phi_1(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\mu},\qquad \phi_2(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w},$$

$$\phi_3(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Phi}(\mathbf{w}\otimes\mathbf{w}),\qquad \phi_4(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Psi}(\mathbf{w}\otimes\mathbf{w}\otimes\mathbf{w}),$$

where $\boldsymbol{\Phi}\in\mathbb{R}^{N\times N^2}$ and $\boldsymbol{\Psi}\in\mathbb{R}^{N\times N^3}$ are flattened co-skewness and co-kurtosis objects. In this convention, skewness and kurtosis are raw third and fourth central moments, not standardized moments (Zhou et al., 2020).

In moment-feasibility and distributional-shape studies, the standard objects are instead the centralized and standardized moments

$D_n=\frac{m_n}{m_2^{n/2}},$

so that $D_3$ is skewness and $D_4$ is kurtosis. A related inferential convention writes

$$\tau=\frac{E[(X-\mu)^3]}{\sigma^3},\qquad \kappa=\frac{E[(X-\mu)^4]}{\sigma^4},\qquad \Delta=\kappa-\tau^2,$$

with $\kappa_e=\kappa-3$ reserved for excess kurtosis. These conventions are location-scale invariant and are therefore natural for feasibility regions in the $$\phi_1(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\mu},\qquad \phi_2(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w},$$0-plane (Meer et al., 2023, Klaassen et al., 2023).

A third convention arises in lattice QCD, where the first four baryon-number cumulants are operationalized through cumulant ratios rather than raw moments. The mean, variance, skewness, and kurtosis of net baryon-number fluctuations are encoded as

$$\phi_1(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\mu},\qquad \phi_2(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w},$$1

$$\phi_1(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\mu},\qquad \phi_2(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w},$$2

Accordingly, the same four-moment problem can be phrased either as an optimization over portfolio weights, a feasibility question for standardized moments, or a cumulant-ratio analysis along constrained thermodynamic trajectories (Collaboration et al., 2017).

2. Optimization formulations in portfolio selection

The central portfolio formulation is a weighted-sum scalarization,

$$\phi_1(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\mu},\qquad \phi_2(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w},$$3

with signs chosen to maximize odd moments and minimize even moments. One widely studied feasible set is

$$\phi_1(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\mu},\qquad \phi_2(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w},$$4

while another standard choice is the simplex

$$\phi_1(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\mu},\qquad \phi_2(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w},$$5

In both cases the resulting problem is a quartic nonconvex polynomial optimization problem. A recurring calibration device chooses $$\phi_1(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\mu},\qquad \phi_2(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w},$$6 from the fourth-order Taylor approximation of CRRA utility, for example

$$\phi_1(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\mu},\qquad \phi_2(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w},$$7

or, equivalently in alternative notation, $$\phi_1(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\mu},\qquad \phi_2(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w},$$8 with the same preference interpretation (Zhou et al., 2020, Niu et al., 2019, Wang et al., 2022).

A distinct but related formulation is the MVSK tilting problem. Rather than minimizing a weighted sum, it seeks the largest common improvement factor $$\phi_1(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\mu},\qquad \phi_2(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w},$$9 relative to a reference portfolio $$\phi_3(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Phi}(\mathbf{w}\otimes\mathbf{w}),\qquad \phi_4(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Psi}(\mathbf{w}\otimes\mathbf{w}\otimes\mathbf{w}),$$0, subject to simultaneous improvement in mean and skewness, simultaneous reduction in variance and kurtosis, a tracking-error bound, and the same portfolio-feasibility constraints. This is a constraint-based improvement model rather than a scalarized utility model (Zhou et al., 2020).

The multi-objective perspective makes the same structure explicit. One may write the four objectives as

$$\phi_3(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Phi}(\mathbf{w}\otimes\mathbf{w}),\qquad \phi_4(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Psi}(\mathbf{w}\otimes\mathbf{w}\otimes\mathbf{w}),$$1

in investor-preference language, or equivalently as minimization of

$$\phi_3(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Phi}(\mathbf{w}\otimes\mathbf{w}),\qquad \phi_4(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Psi}(\mathbf{w}\otimes\mathbf{w}\otimes\mathbf{w}),$$2

Linear scalarization by $$\phi_3(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Phi}(\mathbf{w}\otimes\mathbf{w}),\qquad \phi_4(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Psi}(\mathbf{w}\otimes\mathbf{w}\otimes\mathbf{w}),$$3 produces a quartic polynomial $$\phi_3(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Phi}(\mathbf{w}\otimes\mathbf{w}),\qquad \phi_4(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Psi}(\mathbf{w}\otimes\mathbf{w}\otimes\mathbf{w}),$$4, and sparse extensions further restrict supports by conditions such as $$\phi_3(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Phi}(\mathbf{w}\otimes\mathbf{w}),\qquad \phi_4(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Psi}(\mathbf{w}\otimes\mathbf{w}\otimes\mathbf{w}),$$5 for forbidden sets $$\phi_3(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Phi}(\mathbf{w}\otimes\mathbf{w}),\qquad \phi_4(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Psi}(\mathbf{w}\otimes\mathbf{w}\otimes\mathbf{w}),$$6. This yields cardinality-type sparsity and correlation-exclusion constraints within the same MVSK framework (Steenkamp, 2023).

The “unrestricted” sample-moment large-scale formulation retains the long-only simplex

$$\phi_3(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Phi}(\mathbf{w}\otimes\mathbf{w}),\qquad \phi_4(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Psi}(\mathbf{w}\otimes\mathbf{w}\otimes\mathbf{w}),$$7

but uses the full sample-moment objective without factor restrictions, parametric return-family restrictions, or structural simplifications of the comoments. In that setting,

$$\phi_3(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Phi}(\mathbf{w}\otimes\mathbf{w}),\qquad \phi_4(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Psi}(\mathbf{w}\otimes\mathbf{w}\otimes\mathbf{w}),$$8

with moments computed directly from the centered return matrix $$\phi_3(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Phi}(\mathbf{w}\otimes\mathbf{w}),\qquad \phi_4(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Psi}(\mathbf{w}\otimes\mathbf{w}\otimes\mathbf{w}),$$9, becomes the reference large-universe MVSK problem (Wang et al., 28 Apr 2026).

3. Nonconvexity, convexity, and algorithmic frameworks

The main mathematical difficulty is that $\boldsymbol{\Phi}\in\mathbb{R}^{N\times N^2}$0 and $\boldsymbol{\Phi}\in\mathbb{R}^{N\times N^2}$1 are nonconvex in the portfolio weights. The objective is a cubic-plus-quartic polynomial, the fourth-moment object $\boldsymbol{\Phi}\in\mathbb{R}^{N\times N^2}$2 has size $\boldsymbol{\Phi}\in\mathbb{R}^{N\times N^2}$3, memory cost is $\boldsymbol{\Phi}\in\mathbb{R}^{N\times N^2}$4, evaluating $\boldsymbol{\Phi}\in\mathbb{R}^{N\times N^2}$5 costs $\boldsymbol{\Phi}\in\mathbb{R}^{N\times N^2}$6, and computing its gradient is also $\boldsymbol{\Phi}\in\mathbb{R}^{N\times N^2}$7. The corresponding quartic polynomial minimization problem over a polytope is explicitly described as NP-hard (Zhou et al., 2020, Niu et al., 2019).

Several algorithmic families address this structure.

Paper	Formulation	Main computational device
(Zhou et al., 2020)	weighted-sum MVSK and MVSK tilting	successive convex approximation
(Niu et al., 2019)	quartic MVSK over the simplex	DCA, Boosted-DCA, DC-SOS
(Steenkamp, 2023)	scalarized MVSK on simplex or cube	convexity regions and FISTA
(Wang et al., 2022)	ghMST-parametric MVSK	robust fixed-point acceleration
(Wang et al., 28 Apr 2026)	unrestricted sample-moment MVSK	Yau’s affine-normal descent

Successive convex approximation exploits explicit gradients and Hessians such as

$\boldsymbol{\Phi}\in\mathbb{R}^{N\times N^2}$8

together with the identities

$\boldsymbol{\Phi}\in\mathbb{R}^{N\times N^2}$9

to build efficient local convex surrogates while avoiding frequent expensive high-order moment evaluations (Zhou et al., 2020).

Difference-of-convex programming treats MVSK as a DC program. The projective DC decomposition uses

$\boldsymbol{\Psi}\in\mathbb{R}^{N\times N^3}$0

where $\boldsymbol{\Psi}\in\mathbb{R}^{N\times N^3}$1 is chosen large enough to make both $\boldsymbol{\Psi}\in\mathbb{R}^{N\times N^3}$2 and $\boldsymbol{\Psi}\in\mathbb{R}^{N\times N^3}$3 convex on $\boldsymbol{\Psi}\in\mathbb{R}^{N\times N^3}$4. A more structure-aware alternative is the DC-SOS decomposition, where the quartic polynomial is written as a difference of two convex sums-of-squares polynomials. DCA linearizes the concave part, while Boosted-DCA adds an Armijo-type line search along the DC descent direction $\boldsymbol{\Psi}\in\mathbb{R}^{N\times N^3}$5 (Niu et al., 2019).

Convex scalarization results isolate regions of hyper-parameter space for which the quartic scalarization is actually convex. The Hessian can be written as

$\boldsymbol{\Psi}\in\mathbb{R}^{N\times N^3}$6

so convexity reduces to nonnegativity of the scalar quadratic

$\boldsymbol{\Psi}\in\mathbb{R}^{N\times N^3}$7

A distribution-independent sufficient condition is

$\boldsymbol{\Psi}\in\mathbb{R}^{N\times N^3}$8

equivalently $\boldsymbol{\Psi}\in\mathbb{R}^{N\times N^3}$9, which guarantees convexity on all of $D_n=\frac{m_n}{m_2^{n/2}},$0 (Steenkamp, 2023).

Parametric reformulations replace explicit co-skewness and co-kurtosis tensors by a generalized hyperbolic multivariate skew-$D_n=\frac{m_n}{m_2^{n/2}},$1 model. Under that model,

$D_n=\frac{m_n}{m_2^{n/2}},$2

$D_n=\frac{m_n}{m_2^{n/2}},$3

so memory drops from $D_n=\frac{m_n}{m_2^{n/2}},$4 to $D_n=\frac{m_n}{m_2^{n/2}},$5. The resulting fixed-point formulation uses the projected map

$D_n=\frac{m_n}{m_2^{n/2}},$6

and solves for stationary points through robust fixed-point acceleration (Wang et al., 2022).

Large-scale tensor-free methods work directly with the return matrix $D_n=\frac{m_n}{m_2^{n/2}},$7. In that formulation,

$D_n=\frac{m_n}{m_2^{n/2}},$8

the Hessian factorizes as

$D_n=\frac{m_n}{m_2^{n/2}},$9

and exact line search along any feasible direction reduces to minimizing a scalar quartic. This separates data-map geometry from investor preference coefficients and makes full-universe MVSK feasible in the hundreds and thousands of assets (Wang et al., 28 Apr 2026).

4. Feasibility regions for skewness and kurtosis

A recurrent misconception is that skewness and kurtosis can be chosen freely once mean and variance are fixed. Several papers show that this is false. Pearson’s inequality gives

$D_3$0

equivalently

$D_3$1

For unrestricted distributions with finite fourth moment, this lower bound is exact: the attainable skewness-kurtosis set is

$D_3$2

Under unimodality the restriction strengthens to

$D_3$3

and, more sharply, to a boundary $D_3$4. For symmetric unimodal distributions the exact set becomes

$D_3$5

Thus admissibility of higher moments depends materially on shape restrictions, not only on fourth-moment finiteness (Klaassen et al., 2023).

Support information sharpens these restrictions further. For nonnegative support,

$D_3$6

and for compact support $D_3$7,

$D_3$8

Kurtosis is then subject not only to Pearson’s inequality but also to support-aware lower bounds, including the one-sided formula

$D_3$9

when $D_4$0. These are necessary feasibility constraints for any moment-constrained optimization model (Meer et al., 2023).

An analogous phenomenon appears in concentration theory for bounded martingale differences. If $D_4$1, $D_4$2, and one has conditional skewness lower bounds $D_4$3 or conditional kurtosis upper bounds $D_4$4, then the conditional variance is forced down through

$D_4$5

and

$D_4$6

The resulting effective variance parameter can replace the variance-only term in Hoeffding/Bentkus-type tail bounds (Bentkus et al., 2011).

5. Estimation, inference, and tail-risk translation

Dynamic MVSK analysis requires estimation of conditional higher moments. The quantiled conditional moments method constructs a linear regression across many estimated conditional quantiles,

$D_4$7

with $D_4$8. Regressing estimated conditional quantiles on $D_4$9, $$\tau=\frac{E[(X-\mu)^3]}{\sigma^3},\qquad \kappa=\frac{E[(X-\mu)^4]}{\sigma^4},\qquad \Delta=\kappa-\tau^2,$$0, and $$\tau=\frac{E[(X-\mu)^3]}{\sigma^3},\qquad \kappa=\frac{E[(X-\mu)^4]}{\sigma^4},\qquad \Delta=\kappa-\tau^2,$$1 yields OLS estimates of

$$\tau=\frac{E[(X-\mu)^3]}{\sigma^3},\qquad \kappa=\frac{E[(X-\mu)^4]}{\sigma^4},\qquad \Delta=\kappa-\tau^2,$$2

and hence

$$\tau=\frac{E[(X-\mu)^3]}{\sigma^3},\qquad \kappa=\frac{E[(X-\mu)^4]}{\sigma^4},\qquad \Delta=\kappa-\tau^2,$$3

Under the paper’s assumptions these estimators are consistent and have $$\tau=\frac{E[(X-\mu)^3]}{\sigma^3},\qquad \kappa=\frac{E[(X-\mu)^4]}{\sigma^4},\qquad \Delta=\kappa-\tau^2,$$4 convergence rate (Zhang et al., 2023).

Inference on the skewness-kurtosis pair can also be organized through the scale-free statistic

$$\tau=\frac{E[(X-\mu)^3]}{\sigma^3},\qquad \kappa=\frac{E[(X-\mu)^4]}{\sigma^4},\qquad \Delta=\kappa-\tau^2,$$5

Its empirical estimator is

$$\tau=\frac{E[(X-\mu)^3]}{\sigma^3},\qquad \kappa=\frac{E[(X-\mu)^4]}{\sigma^4},\qquad \Delta=\kappa-\tau^2,$$6

and, under finite eighth moment,

$$\tau=\frac{E[(X-\mu)^3]}{\sigma^3},\qquad \kappa=\frac{E[(X-\mu)^4]}{\sigma^4},\qquad \Delta=\kappa-\tau^2,$$7

This supports tests of the form $$\tau=\frac{E[(X-\mu)^3]}{\sigma^3},\qquad \kappa=\frac{E[(X-\mu)^4]}{\sigma^4},\qquad \Delta=\kappa-\tau^2,$$8, which detect non-unimodality when the estimated skewness-kurtosis pair falls below the unimodal lower bound (Klaassen et al., 2023).

The first four moments can also be translated into downside-risk surrogates. For a loss variable $$\tau=\frac{E[(X-\mu)^3]}{\sigma^3},\qquad \kappa=\frac{E[(X-\mu)^4]}{\sigma^4},\qquad \Delta=\kappa-\tau^2,$$9 with positive skewness $\kappa_e=\kappa-3$0 and positive excess kurtosis $\kappa_e=\kappa-3$1, one proposed VaR approximation is

$\kappa_e=\kappa-3$2

with a matching ES approximation defined by an integral correction to the Gaussian ES term. This provides a direct mapping from $\kappa_e=\kappa-3$3 to tail-risk proxies, although the paper also states that heavy-tail cases such as Pareto and lognormal can produce substantial approximation error as $\kappa_e=\kappa-3$4 (Barczy et al., 2018).

6. Extensions, applications, and limitations

Although portfolio selection is the main applied setting, the same four-moment structure appears in lattice QCD. There the focus is not optimization over weights but Taylor expansion of cumulant ratios at small baryon chemical potential under strangeness neutrality and fixed charge-to-baryon ratio. At fixed temperature,

$\kappa_e=\kappa-3$5

$\kappa_e=\kappa-3$6

so the mean ratio is odd in $\kappa_e=\kappa-3$7, while the skewness and kurtosis ratios are even. This suggests that the mean-variance-skewness-kurtosis problem is also a structured cumulant-comparison problem in statistical physics, not only a financial optimization problem (Collaboration et al., 2017).

Recent large-scale portfolio evidence indicates that unrestricted sample-moment MVSK can be solved directly on very large universes and compared against exact mean-variance portfolios. On a 5-minute A-share panel with 5,440 stocks, the reported out-of-sample comparisons show that the incremental value of higher moments is strongest at moderate return targets, while the downside advantage largely disappears at the most aggressive targets. This suggests that the practical value of skewness and kurtosis control depends materially on the return target and the residual design freedom left after the mean constraint is imposed (Wang et al., 28 Apr 2026).

Several limitations are recurrent across the literature. Weighted-sum scalarization recovers only part of the Pareto front in general, even when many scalarizations are solved on a dense grid. Feasibility results show that not every target $\kappa_e=\kappa-3$8 or $\kappa_e=\kappa-3$9 is admissible. The unimodality criterion based on $$\phi_1(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\mu},\qquad \phi_2(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w},$$00 is one-sided: large $$\phi_1(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\mu},\qquad \phi_2(\mathbf{w})=\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w},$$01 does not imply unimodality. Four-moment VaR and ES approximations can be informative, but they remain approximations and may deteriorate for strongly non-Gaussian tails. Accordingly, the mean-variance-skewness-kurtosis problem is best understood not as a single optimization template, but as a family of higher-moment models whose validity depends on objective choice, distributional assumptions, admissibility constraints, and computational regime (Steenkamp, 2023, Klaassen et al., 2023, Barczy et al., 2018).