# Expressions for moments of sample cumulants?

I sample $$n$$ points from standard normal and need the mean and variance of 3rd and 4th sample cumulants.

@JimB suggested that variance of 4th sample cumulant is given by the expression below. This seems like a very algebraically nice expression, is there way to derive/verify it with the help of Mathematica?

$$\frac {24(n-1)((n-6)n+24)}{n^4}$$

Edit, Sep 5

The following gives a way to verify the formula for a fixed $$n$$, but not helpful for deriving the formula in terms of $$n$$

n = 5;
Clear[v];
rvec = Array[v, n];
sampleDist = ProductDistribution @@ Table[NormalDistribution[], n];
cumulantDist =
TransformedDistribution[Cumulant[rvec, 4],
rvec \[Distributed] sampleDist];
Variance[cumulantDist] == (24 (n - 1) ((n - 6) n + 24))/n^4

• Have a look at the mathstatica package, they also have a free book mathstatica.com/book/bookcontents.html Sep 6, 2020 at 21:03
• @yarchik thanks for the pointer, it seems this issue is discussed in Section 7.3 "moments of moments". I actually bought 1.0 version of their package a while back, good to see they are still maintaining it Sep 7, 2020 at 4:32

Here is a solution which is completely general, and works for any random variable whose moments exist. The standard Normal is then just a special and simple case.

The OP is interested in the $$3^\text{rd}$$ and $$4^\text{th}$$ sample cumulants (a type of sample moment, similar to sample central moments). The lingua franca for working with sample moments, and moments of sample moments, is to structure everything in terms of power sum notation, $$s_r = \sum_{i=1}^n X_i^r$$.

The mathStatica package for Mathematica has a function to express sample cumulants into power sum notation. Here is the $$4^\text{th}$$ sample cumulant, expressed in terms of power sums: This statistic is a random variable, and the OP wishes to find some moments of it.

The mean of k4 is simply the $$1^\text{st}$$ Raw Moment of k4: The variance of k4 is the $$2^\text{nd}$$ Central Moment of k4: ... where the solution is expressed in terms of the central moments $$\mu_i$$ of the population of the parent random variable $$X$$.

Note that this solution is completely general, and holds for any random variable whose moments exist.

## The standard Normal case

In the simple standard Normal case $$X \sim N(0,1)$$, the first 8 central moments $$\mu_r$$ of the population are: Substituting these values into the solutions above yields the mean of k4 as: and the variance of k4 as: All done.

This is similar approach using FindSequenceFunction using the sample cumulant definition that Mathematica uses. (That sample cumulant definition is not always unbiased for the population cumulant.)

First define a function that finds the variance for the k-th cumulant for a sample size of $$n$$ from a Normal[0,1] distribution.

variance[k_, n_] := Module[{e1, e2, expectations},
(* Finds variance of k-th cumulant for a sample of size n from a Normal[0,1] distribution *)
(* Raw moments of a Normal[0,1] *)
expectations = Table[x[i_]^i -> Expectation[z^i, z \[Distributed] NormalDistribution[0, 1]],
{i, 2 k, 1, -1}];
(* Expectation of sample cumulant *)
e1 = (Cumulant[Table[x[i], {i, n}], k] // Expand) /. expectations;
(* Expecation of square of sample cumlant *)
e2 = (Cumulant[Table[x[i], {i, n}], k]^2 // Expand) /. expectations;
e2 - e1^2]


Now find the variance for the 3rd and 4th cumulants:

v3 = Table[{n, variance[3, n]}, {n, 2, 10}];
FindSequenceFunction[v3, n] // FullSimplify


$$\frac{6 (n-2) (n-1)}{n^3}$$

v4 = Table[{n, variance[4, n]}, {n, 2, 10}];
FindSequenceFunction[v4, n] // FullSimplify


$$\frac{24 (n-1) ((n-6) n+24)}{n^4}$$

mathStatica (for me) has the more straightforward way to getting what you asked (as @wolfies has demonstrated). Here is how to accomplish that in a somewhat similar fashion in Mathematica.

(* Sample statistic *)
(k4 = MomentConvert[Cumulant, "SampleEstimator"] /. PowerSymmetricPolynomial -> n) // TraditionalForm (* Mean of sample statistic *)
(meank4 = MomentConvert[k4, "CentralMoment"] /. PowerSymmetricPolynomial -> n // FullSimplify) // TraditionalForm (* Variance of sample statistic *)
(vark4 = MomentConvert[k4^2 - meank4^2, "CentralMoment"] /. PowerSymmetricPolynomial -> n // Simplify) // TraditionalForm These can be simplified under the assumption of random samples from a Normal(0,1) distribution:

meank4 /. CentralMoment[i_] -> CentralMoment[NormalDistribution[0, 1], i] // FullSimplify vark4 /. CentralMoment[i_] -> CentralMoment[NormalDistribution[0, 1], i] // FullSimplify An idea from JimB was to use FindSequenceFunction

getExpr[n_] := Module[{v},
rvec = Array[v, n];
sampleDist = ProductDistribution @@ Table[NormalDistribution[], n];
cumulantDist =
TransformedDistribution[Cumulant[rvec, 4],
rvec \[Distributed] sampleDist]; Variance[cumulantDist]
];
values = Table[{n, getExpr[n]}, {n, 1, 9}];
FindSequenceFunction[values] (* (24 (-24 + 30 #1 - 7 #1^2 + #1^3))/#1^4 & *)


An alternative way to use MomentConvert + MomentEvaluate:

## $$\kappa _3$$

Mean

μκ3 = FullSimplify @
MomentConvert[Cumulant, {"SampleEstimator", n}, {"CentralMoment", n}]; MomentEvaluate[μκ3, NormalDistribution[]]

0


Variance

σκ3 = FullSimplify[
MomentConvert[Cumulant^2, {"SampleEstimator", n}, {"CentralMoment", n}] - μκ3^2]; MomentEvaluate[σκ3, NormalDistribution[]] ## $$\kappa _4$$

Mean

μκ4 = FullSimplify @
MomentConvert[Cumulant, {"SampleEstimator", n}, {"CentralMoment", n}]; MomentEvaluate[μκ4, NormalDistribution[]] Variance

σκ4 = FullSimplify[
MomentConvert[Cumulant^2, {"SampleEstimator", n}, {"CentralMoment", n}] - μκ4^2]; MomentEvaluate[σκ4, NormalDistribution[]] 