Covariance of sample central moments

Question

What is the covariance matrix of the first four sample central moments?
Assuming a zero mean and possibly without assuming normality.

The covariance matrix of the first four sample raw moments for a normal random variable of zero mean is:

$$ \left( \begin{array}{cccc} s^2 & 0 & 3 s^4 & 0 \\ 0 & 2 s^4 & 0 & 12 s^6 \\ 3 s^4 & 0 & 15 s^6 & 0 \\ 0 & 12 s^6 & 0 & 96 s^8 \\ \end{array} \right) $$ and can be obtained e.g. in Mathematica using the command
Covariance[ TransformedDistribution[{x, x^2, x^3, x^4}, Distributed[x, NormalDistribution[0, s]]]]

The variance of the sample variance was discussed here.

Answer 1

OP asks: What is the covariance matrix of the first four sample central moments?

The Question: more formally

Let $(X_1, ..., X_n)$ denote a random sample of size $n$ drawn from a population random variable $X$. Then, by definition, the $r$th sample central moment is:

$$m_r = \frac1n \sum _{i=1}^n \left(X_i-\bar{X}\right){}^r$$

The problem is to find: $\operatorname{cov}(m_i, m_j)$ for $i,j$ = 1 to 4.

Background

The problem is just a special case of the much more general problem of 'moments of moments' which are usually defined in terms of power sum notation:

$$s_r=\sum _{i=1}^n X_i^r$$

There are those on math.SE who like to see solutions posted with every step shown ... old-fashioned long-division-style. The question posed here is somewhat ill-suited to step-by-step derivation by hand, so much so, that doing it by hand has resulted in the published literature in the field being riddled with errors and mistakes. Indeed, there are even some rather famous papers such as Fisher's (1928): "Moments and product moments of sampling distributions", Proceedings of the London Mathematical Society - series 2, volume 30 ... the results of which were subsequently published across the 6 editions of reference bibles such as Kendall's Advanced Theory of Statistics ... that for more than 70 years remained significantly in error until they were corrected by the mathStatica software in 2002 (of which I am an author).

Fortunately today, with the aid of automated computer algebra, these problems are now relatively straightforwards, leaving the grunt work to your computer ...

Solution

There are basically two steps:

Step 1

The first step is to express the sample central moments $m_r$ in terms of power sums $s_i$. This is done here by the mathStatica function SampleCentralToPowerSum[r]:

Step 2

Since $\operatorname{cov}(m_i, m_j)$ = $\mu_{1,1}(m_i, m_j)$ ... i.e. since the covariance operator is just the {1,1} CentralMoment ... the solution can be obtained with the mathStatica function:

CentralMomentToCentral[{1,1}, {m_i, m_j}]

One could now pop out the whole 4x4 covariance matrix desired as a one-liner, but the solutions are quite complicated expressions, and the full matrix will not fit on screen, so it is better here to present them row by row ...

• First Row: since $m_1 = 0$, $\operatorname{cov}(m_1, m_j)$ = 0

• Second Row:

where $\mu_i$ denotes the $i$th central moment of the population of $X$.

• Third Row:

• Fourth Row:

All other results by symmetry.

All done.

---

Other: Comment on OP's 'solution' using Mathematica

The OP writes:

The covariance matrix of the first four central moments for a normal random variable is:

<p>$$

\left( \begin{array}{cccc} s^2 & 0 & 3 s^4 & 0 \ 0 & 2 s^4 & 0 & 12 s^6 \ 3 s^4 & 0 & 15 s^6 & 0 \ 0 & 12 s^6 & 0 & 96 s^8 \ \end{array} \right) $$ and can be obtained e.g. in Mathematica using the command
Covariance[ TransformedDistribution[{x, x^2, x^3, x^4}, Distributed[x, NormalDistribution[0, s]]]]

---

This is incorrect for a variety of reasons:

• First, you assert that you are calculating the covariance matrix for the first 4 CENTRAL moments. Please note that the first 4 Central moments are constants ... not random variables ... so the covariance matrix of the first 4 central moments is a null matrix ... which is obviously not what you want.

• Second, there is nothing central about your calculation in Mathematica, since you have set the mean mu to 0, ... so it is all raw.

• Third, what you are actually calculating here is $\operatorname{cov}(X^i, X^j)$ for $i,j$ = 1 to 4 when $X$ ~ $N(0,s^2)$. This is not the problem posed at the top, namely to find $\operatorname{cov}(m_i, m_j)$ for $i,j$ = 1 to 4.

Answer 2

The following is from jupyter notebook by python:

</p>

In [1]:

from sympy import S, sqrt
from SampleStats import SampleStats
from utils import to_latex
from iprint import *

# for auto-reloading extenrnal modules
%load_ext autoreload
%autoreload 2

In [2]:

ss = SampleStats()
n = ss.n

In [3]:

iprint_eqnarray([(k, v) for k, v in ss.moments_c2r_map.items()][:8])

$\displaystyle \begin{align} \mu_{2} &= - \mu^{2} + \mu'{2}, \ \mu{3} &= 2 \mu^{3} - 3 \mu \mu'{2} + \mu'{3}, \ \mu_{4} &= - 3 \mu^{4} + 6 \mu^{2} \mu'{2} - 4 \mu \mu'{3} + \mu'{4}, \ \mu{5} &= 4 \mu^{5} - 10 \mu^{3} \mu'{2} + 10 \mu^{2} \mu'{3} - 5 \mu \mu'{4} + \mu'{5}, \ \mu_{6} &= - 5 \mu^{6} + 15 \mu^{4} \mu'{2} - 20 \mu^{3} \mu'{3} + 15 \mu^{2} \mu'{4} - 6 \mu \mu'{5} + \mu'{6}, \ \mu{7} &= 6 \mu^{7} - 21 \mu^{5} \mu'{2} + 35 \mu^{4} \mu'{3} - 35 \mu^{3} \mu'{4} + 21 \mu^{2} \mu'{5} - 7 \mu \mu'{6} + \mu'{7}, \ \mu_{8} &= - 7 \mu^{8} + 28 \mu^{6} \mu'{2} - 56 \mu^{5} \mu'{3} + 70 \mu^{4} \mu'{4} - 56 \mu^{3} \mu'{5} + 28 \mu^{2} \mu'{6} - 8 \mu \mu'{7} + \mu'{8}, \ \mu{9} &= 8 \mu^{9} - 36 \mu^{7} \mu'{2} + 84 \mu^{6} \mu'{3} - 126 \mu^{5} \mu'{4} + 126 \mu^{4} \mu'{5} - 84 \mu^{3} \mu'{6} + 36 \mu^{2} \mu'{7} - 9 \mu \mu'{8} + \mu'{9}. \end{align}$

In [4]:

iprint_eqnarray([(k, v) for k, v in ss.moments_r2c_map.items()][:8])

$\displaystyle \begin{align} \mu'{2} &= \mu^{2} + \mu{2}, \ \mu'{3} &= \mu^{3} + 3 \mu \mu{2} + \mu_{3}, \ \mu'{4} &= \mu^{4} + 6 \mu^{2} \mu{2} + 4 \mu \mu_{3} + \mu_{4}, \ \mu'{5} &= \mu^{5} + 10 \mu^{3} \mu{2} + 10 \mu^{2} \mu_{3} + 5 \mu \mu_{4} + \mu_{5}, \ \mu'{6} &= \mu^{6} + 15 \mu^{4} \mu{2} + 20 \mu^{3} \mu_{3} + 15 \mu^{2} \mu_{4} + 6 \mu \mu_{5} + \mu_{6}, \ \mu'{7} &= \mu^{7} + 21 \mu^{5} \mu{2} + 35 \mu^{4} \mu_{3} + 35 \mu^{3} \mu_{4} + 21 \mu^{2} \mu_{5} + 7 \mu \mu_{6} + \mu_{7}, \ \mu'{8} &= \mu^{8} + 28 \mu^{6} \mu{2} + 56 \mu^{5} \mu_{3} + 70 \mu^{4} \mu_{4} + 56 \mu^{3} \mu_{5} + 28 \mu^{2} \mu_{6} + 8 \mu \mu_{7} + \mu_{8}, \ \mu'{9} &= \mu^{9} + 36 \mu^{7} \mu{2} + 84 \mu^{6} \mu_{3} + 126 \mu^{5} \mu_{4} + 126 \mu^{4} \mu_{5} + 84 \mu^{3} \mu_{6} + 36 \mu^{2} \mu_{7} + 9 \mu \mu_{8} + \mu_{9}. \end{align}$

In [5]:

p = 4
eqn_list = []
for i in range(1, p+1):
    iprint_eqs('\mathbb{{E}}[\\bar{{X}}^{{{:d}}}]'.format(i),
               to_latex(ss.E(ss.mean**i)))
    print()

$\displaystyle \mathbb{E}[\bar{X}^{1}]=\mu$

$\displaystyle \mathbb{E}[\bar{X}^{2}]=\mu^{2} + \frac{\mu_{2}}{n}$

$\displaystyle \mathbb{E}[\bar{X}^{3}]=\mu^{3} + \frac{3 \mu \mu_{2}}{n} + \frac{\mu_{3}}{n^{2}}$

$\displaystyle \mathbb{E}[\bar{X}^{4}]=\mu^{4} + \frac{6 \mu^{2} \mu_{2}}{n} + \frac{4 \mu \mu_{3}}{n^{2}} + \frac{3 \mu_{2}^{2} \left(n - 1\right)}{n^{3}} + \frac{\mu_{4}}{n^{3}}$

In [6]:

M = [ss.central_moment(i) for i in range(p+1)]
eqn_list = []
for i in range(1, p+1):
    iprint_eqs('\mathbb{{E}}[M_{{{:d}}}]'.format(i),
               to_latex(M[i]), to_latex(ss.E(M[i])))
    print()

$\displaystyle \begin{align} \mathbb{E}[M_{1}] &= 0 \ &= 0. \end{align}$

$\displaystyle \begin{align} \mathbb{E}[M_{2}] &= - \frac{S_{1}^{2}}{n^{2}} + \frac{S_{2}}{n} \ &= \frac{\mu_{2} \left(n - 1\right)}{n}. \end{align}$

$\displaystyle \begin{align} \mathbb{E}[M_{3}] &= \frac{2 S_{1}^{3}}{n^{3}} - \frac{3 S_{1} S_{2}}{n^{2}} + \frac{S_{3}}{n} \ &= \frac{\mu_{3} \left(n - 2\right) \left(n - 1\right)}{n^{2}}. \end{align}$

$\displaystyle \begin{align} \mathbb{E}[M_{4}] &= - \frac{3 S_{1}^{4}}{n^{4}} + \frac{6 S_{1}^{2} S_{2}}{n^{3}} - \frac{4 S_{1} S_{3}}{n^{2}} + \frac{S_{4}}{n} \ &= \frac{3 \mu_{2}^{2} \left(n - 1\right) \left(2 n - 3\right)}{n^{3}} + \frac{\mu_{4} \left(n - 1\right) \left(n^{2} - 3 n + 3\right)}{n^{3}}. \end{align}$

In [7]:

iprint_eqnarray([('\mathrm{{Var}}[M_{{{:d}}}]'.format(i),
                  to_latex(ss.var(M[i])))
                 for i in range(1, p+1)])

$\displaystyle \begin{align} \mathrm{Var}[M_{1}] &= 0, \ \mathrm{Var}[M_{2}] &= - \frac{\mu_{2}^{2} \left(n - 3\right) \left(n - 1\right)}{n^{3}} + \frac{\mu_{4} \left(n - 1\right)^{2}}{n^{3}}, \ \mathrm{Var}[M_{3}] &= \frac{3 \mu_{2}^{3} \left(n - 2\right) \left(n - 1\right) \left(3 n^{2} - 12 n + 20\right)}{n^{5}} - \frac{3 \mu_{2} \mu_{4} \left(n - 2\right)^{2} \left(n - 1\right) \left(2 n - 5\right)}{n^{5}} - \frac{\mu_{3}^{2} \left(n - 10\right) \left(n - 2\right)^{2} \left(n - 1\right)}{n^{5}} + \frac{\mu_{6} \left(n - 2\right)^{2} \left(n - 1\right)^{2}}{n^{5}}, \ \mathrm{Var}[M_{4}] &= \frac{18 \mu_{2}^{4} \left(n - 1\right) \left(4 n^{4} - 52 n^{3} + 228 n^{2} - 438 n + 315\right)}{n^{7}} - \frac{12 \mu_{2}^{2} \mu_{4} \left(n - 1\right) \left(10 n^{4} - 93 n^{3} + 324 n^{2} - 513 n + 315\right)}{n^{7}} + \frac{8 \mu_{2} \mu_{3}^{2} \left(n - 2\right) \left(n - 1\right) \left(2 n^{4} - 24 n^{3} + 120 n^{2} - 315 n + 315\right)}{n^{7}} + \frac{4 \mu_{2} \mu_{6} \left(n - 1\right) \left(n^{2} - 3 n + 3\right) \left(10 n^{2} - 27 n + 21\right)}{n^{7}} - \frac{8 \mu_{3} \mu_{5} \left(n - 1\right) \left(n^{2} - 3 n + 3\right) \left(n^{3} - 4 n^{2} + 18 n - 21\right)}{n^{7}} - \frac{\mu_{4}^{2} \left(n - 1\right) \left(n^{5} - 23 n^{4} + 117 n^{3} - 345 n^{2} + 531 n - 315\right)}{n^{7}} + \frac{\mu_{8} \left(n - 1\right)^{2} \left(n^{2} - 3 n + 3\right)^{2}}{n^{7}}. \end{align}$

In [8]:

iprint_eqnarray([('\mathrm{{Cov}}[M_{{{:d}}}, M_{{{:d}}}]'.format(i, j),
                  to_latex(ss.cov(M[i], M[j])))
                 for i in range(1, p+1)
                 for j in range(1, p+1) if i < j])

$\displaystyle \begin{align} \mathrm{Cov}[M_{1}, M_{2}] &= 0, \ \mathrm{Cov}[M_{1}, M_{3}] &= 0, \ \mathrm{Cov}[M_{1}, M_{4}] &= 0, \ \mathrm{Cov}[M_{2}, M_{3}] &= - \frac{2 \mu_{2} \mu_{3} \left(n - 2\right) \left(n - 1\right) \left(2 n - 5\right)}{n^{4}} + \frac{\mu_{5} \left(n - 2\right) \left(n - 1\right)^{2}}{n^{4}}, \ \mathrm{Cov}[M_{2}, M_{4}] &= - \frac{6 \mu_{2}^{3} \left(n - 1\right) \left(4 n^{2} - 15 n + 15\right)}{n^{5}} - \frac{\mu_{2} \mu_{4} \left(n - 1\right) \left(n^{3} - 24 n^{2} + 60 n - 45\right)}{n^{5}} - \frac{2 \mu_{3}^{2} \left(n - 1\right) \left(2 n^{3} - 8 n^{2} + 18 n - 15\right)}{n^{5}} + \frac{\mu_{6} \left(n - 1\right)^{2} \left(n^{2} - 3 n + 3\right)}{n^{5}}, \ \mathrm{Cov}[M_{3}, M_{4}] &= \frac{6 \mu_{2}^{2} \mu_{3} \left(n - 2\right) \left(n - 1\right) \left(2 n^{3} - 24 n^{2} + 86 n - 105\right)}{n^{6}} - \frac{3 \mu_{2} \mu_{5} \left(n - 2\right) \left(n - 1\right) \left(n^{3} - 12 n^{2} + 28 n - 21\right)}{n^{6}} - \frac{\mu_{3} \mu_{4} \left(n - 2\right) \left(n - 1\right) \left(5 n^{3} - 32 n^{2} + 102 n - 105\right)}{n^{6}} + \frac{\mu_{7} \left(n - 2\right) \left(n - 1\right)^{2} \left(n^{2} - 3 n + 3\right)}{n^{6}}. \end{align}$

In [9]:

G1 = sqrt(n * (n-1)) / (n-2) * ss.E(M[3]) / ss.E(M[2])**(S(3)/2)
G1.simplify()

Out[9]:

$\displaystyle \frac{\mu_{3}}{\mu_{2}^{\frac{3}{2}}}$

In [10]:

G2 = ss.E((n-1) * ((n+1)*M[4] - 3*(n-1)*M[2]**2)) / ((n-2) * (n-3) * ss.E(M[2])**2)
G2.simplify()

Out[10]:

$\displaystyle -3 + \frac{\mu_{4}}{\mu_{2}^{2}}$

In [ ]:

</div>