Given a sample of a real-valued time series, $x=\{x_t\}_{t=1,...,T}$, let $\bar{x}$ be the sample mean of $x$ and set $y_t=(x_t-\bar{x})^2$. Then, $\bar{y}$ estimates the variance of $x$.

Question: What is the covariance of $\bar{x}$ and $\bar{y}$? Can I simply take $\frac{1}{T}\sum\limits_{t=1}^T (x_t-\bar{x})(y_t-\bar{y})$?

[Note: I ignore the $\frac{1}{T-1}$ bias correction as it does not arise in MLE or GMM applications.]


2 Answers 2


The OP seeks the covariance (i.e. the $\mu_{1,1}$ central moment) of two sample moments. The question posed is a relatively simple example of a much more general problem known as finding moments of moments. While these problems can rapidly involve lots of messy algebra (even in this simple case), the use of computer algebra systems can make the process easy. The modus operandi for solving such problems is to work with power sum notation $s_r$, namely:

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

In this case, the OP is interested in the sample mean $\frac{s_1}{n}$, and the $2^\text{nd}$ sample central moment expressed in power sums as $m_2 = \frac{1}{n}\sum_{i=1}^n \left(x_i - \bar x\right)^2 = \frac{1}{n}\left(s_2 - \frac{s^2_1}{n}\right).$

Solution (one-liner)

The covariance operator is just the $\mu_{1,1}$ central moment ... so $\text{Cov}(\frac{s_1}{n},m_2)$ can be found using the CentralMomentToCentral function from the mathStatica package for Mathematica as:

enter image description here

where $\mu_r$ denotes the $r^{th}$ central moment of $X$ i.e. $\mu_3 = E[(X-\mu)^3]$. The answer is the same as that given by Ben (after simplifying out skewness), and holds for any distribution whose moments exist.

General solution for 2 variables: variance-covariance matrix

In the comments, the OP asks the more general question:

Suppose there are two time series ${\{x_t\}_{𝑡=1,...,n}}$ and ${\{y_t\}_{t=1,...,n}}$ with corresponding sample means (for $X$ and $Y$) and sample variances (for $X$ and $Y$). Can we also compute the 4×4 covariance matrix of the 4 sample estimators?

The same tools work just as simply in a bivariate world. Note that we do not need to assume that the variables are independent. Using the same power sum notation, now in a bivariate world, let $s_{r,w}=\sum _{i=1}^n X_i^r Y_i^w$. For example, $s_{1,0} = \sum _{i=1}^n X_i$ and $s_{0,1} =\sum _{i=1}^n Y_i$.

The OP's sample variance estimator for $X$ in terms of bivariate power sums is:

enter image description here

and similarly the OP's sample variance estimator for $Y$ in terms of bivariate power sums is:

enter image description here

The list of our 4 estimators is then:

enter image description here

Then the variance-covariance matrix of the 4 estimators is:

enter image description here


  • $\mu _{r,s}$ denotes the product central moment:

$$\mu _{r,s}=E\left[(X-E[X]]^r (Y-E[Y])^s\right]$$

For example, $\mu_{1,1} = \text{Cov}(X,Y)$, $\mu_{2,0}= \text{Var}(X)$ and $\mu_{0,2}= \text{Var}(Y)$.

This is a general solution for any bivariate distribution whose moments exist.


Firstly, you need to distinguish between the true covariance and the sample estimator for this covariance. Generally we examine these types of cases by first deriving the true covariance, and then seeing if the form of this quantity suggests any natural estimator. We can then examine the properties of various plausible estimators to see which constitute a reasonable sample estimator.

I will show you how to derive the covariance for the case of IID data. For time-series data the IID assumption often does not hold, but we need to walk before we can run. Extension to broader model forms, such as using autocorrelated data, is more onerous, and I leave it as an extension for you if you are interested. A scaled version of this result is given in O'Neill (2014) (Result 8, p. 284), and you can find a range of other useful moment results for these types of quantities in that paper.

Covariance in the IID model: We will assume that we are dealing with the simple case where we have IID observations $X_1,...,X_T$ from a distribution with mean $\mu$, variance $\sigma^2$ and skewness $\gamma$ (we assume that the first three raw moments of the distribution are finite). We begin by writing the quantity $\bar{Y}$ as:

$$\begin{align} \bar{Y} &= \frac{1}{T} \sum_t (X_t - \bar{X})^2 \\[6pt] &= \frac{1}{T} \sum_t (X_t^2 - 2 X_t \bar{X} + \bar{X}^2) \\[6pt] &= \frac{1}{T} \bigg[ \sum_t X_t^2 - \frac{2}{T} \sum_t \sum_r X_t X_r + \frac{1}{T} \sum_s \sum_r X_s X_r \bigg] \\[6pt] &= \frac{1}{T} \bigg[ \sum_t X_t^2 - \frac{1}{T} \sum_t \sum_r X_t X_r \bigg]. \\[6pt] \end{align}$$

This allows us to write the covariance of interest as:

$$\begin{align} \mathbb{Cov}(\bar{X}, \bar{Y}) &= \frac{1}{T^2} \bigg[ \sum_t \sum_i \mathbb{Cov}(X_t^2, X_i) - \frac{1}{T} \sum_t \sum_r \sum_i \mathbb{Cov}(X_t X_r, X_i) \bigg]. \\[6pt] \end{align}$$

Now, since the values are assumed to be IID, the parts of this expression are given by:

$$\begin{align} \sum_t \sum_i \mathbb{Cov}(X_t^2, X_i) &= \sum_t \mathbb{Cov}(X_t^2, X_t) \\[6pt] &= \sum_t \bigg[ \mathbb{E}(X_t^3) - \mathbb{E}(X_t^2) \mathbb{E}(X_t) \bigg] \\[6pt] &= \sum_t \bigg[ (\gamma \sigma^3 + 3 \mu \sigma^2 + \mu^3) - \mu (\sigma^2 + \mu^2) \bigg] \\[6pt] &= \sum_t \bigg[ \gamma \sigma^3 + 2 \mu \sigma^2 \bigg] \\[10pt] &= T (\gamma \sigma^3 + 2 \mu \sigma^2), \\[16pt] \sum_t \sum_r \sum_i \mathbb{Cov}(X_t X_r, X_i) &= \sum_t \sum_i \mathbb{Cov}(X_t^2, X_i) + \sum_t \sum_{r \neq t} \sum_i \mathbb{Cov}(X_t X_r, X_i) \\[6pt] &= T (\gamma \sigma^3 + 2 \mu \sigma^2) + 2 \sum_t \sum_{r \neq t} \mathbb{Cov}(X_t X_r, X_t) \\[6pt] &= T (\gamma \sigma^3 + 2 \mu \sigma^2) + 2 \sum_t \sum_{r \neq t} \mathbb{E}(X_r) \mathbb{V}(X_t) \\[6pt] &= T (\gamma \sigma^3 + 2 \mu \sigma^2) + 2 \sum_t \sum_{r \neq t} \mu \sigma^2 \\[6pt] &= T (\gamma \sigma^3 + 2 \mu \sigma^2) + 2 T(T-1) \mu \sigma^2 \\[8pt] &= T (\gamma \sigma^3 + 2T \mu \sigma^2). \\[6pt] \end{align}$$

(For a similar moment result, see this related question.) Putting things together then gives the covariance expression:

$$\begin{align} \mathbb{Cov}(\bar{X}, \bar{Y}) &= \frac{1}{T^2} \bigg[ \sum_t \sum_i \mathbb{Cov}(X_t^2, X_i) - \frac{1}{T} \sum_t \sum_r \sum_i \mathbb{Cov}(X_t X_r, X_i) \bigg] \\[6pt] &= \frac{1}{T^2} \bigg[ T (\gamma \sigma^3 + 2 \mu \sigma^2) - (\gamma \sigma^3 + 2T \mu \sigma^2) \bigg] \\[6pt] &= \frac{1}{T^2} \bigg[ T \gamma \sigma^3 + 2 T \mu \sigma^2 - \gamma \sigma^3 - 2T \mu \sigma^2 \bigg] \\[6pt] &= \frac{T-1}{T^2} \cdot \gamma \sigma^3, \\[6pt] \end{align}$$

A scaled version of this result is given in O'Neill (2014) (Result 8, p. 284). As you can see, this expression shows that the covariance of these quantities is largely determined by the skewness of the underlying distribution. In particular, the covariance is zero when the underlying distribution is unskewed. A natural sample estimator for the above covariance quantity is:

$$\hat{\mathbb{Cov}}(\bar{X}, \bar{Y}) = \frac{T-1}{T^2} \cdot K_T \cdot S_T^3,$$

where $K_T$ is an appropriate measure of sample skewness (there are a few of them) and $S_T^2$ is the sample variance. This will give you a consistent estimator for the true covariance, and the estimator should perform reasonably well in most cases.

Correlation coefficient: You might also be interested to know the corresponding correlation coefficient in this situation. Applying the known variances of the sample mean and sample variance (see O'Neill (2014), Results 2-3, p. 284) we get:

$$\begin{align} \mathbb{Corr}(\bar{X}, \bar{Y}) &= \frac{\mathbb{Cov}(\bar{X}, \bar{Y})}{\mathbb{S}(\bar{X}) \cdot \mathbb{S}(\bar{Y})} \\[6pt] &= \frac{T-1}{T^2} \cdot \gamma \sigma^3 \Bigg/ \frac{\sigma}{\sqrt{T}} \cdot \frac{T-1}{T} \sqrt{\kappa - \frac{T-3}{T-1}} \frac{\sigma^2}{\sqrt{T}} \\[6pt] &= \frac{T-1}{T^2} \cdot \gamma \sigma^3 \Bigg/ \frac{T-1}{T^2} \sqrt{\kappa - \frac{T-3}{T-1}} \sigma^3 \\[6pt] &= \frac{\gamma}{\sqrt{\kappa - (T-3)/(T-1)}} . \\[6pt] \end{align}$$

As $T \rightarrow \infty$ we get:

$$\mathbb{Corr}(\bar{X}, \bar{Y}) \rightarrow \frac{\gamma}{\sqrt{\kappa - 1}},$$

which is the adjusted skewness value. This result is also a scaled version of results giving in O'Neill (2014) (Result 12, p. 285). You can see from these results that the direction of the correlation is determined by the skewness. The intuitive reason for this is quite simple --- if the sample variance increases, this will generally occur because there are more extreme values in the tails of the distribution. For positively skewed distribution the positive deviations for these extreme values will tend to outweigh the negative deviations, so the mean will tend to be higher in this case. Contrarily, for negatively skewed distribution the negative deviations for these extreme values will tend to outweigh the positive deviations, so the mean will tend to be lower in this case.

  • $\begingroup$ Thanks a lot. That's a thorough answer. I'm grateful for you providing all the details!:) Does it matter what sample skewness estimator I use for $K_T$? I guess I can deal with non IID using a kernel approach as in Newey West for heteroskedasticity and autocorrelation consistent standard errors. I struggle to come up with an intuitive idea why skewness determines this covariance though. Is there some insight why the estimators are uncorrelated for symmetric distributions? $\endgroup$
    – Alex
    Aug 31, 2021 at 0:37
  • 1
    $\begingroup$ The intuitive reason that skewness matters is that you are seeking a covariance between a first-order value $X$ and a second-order value $Y$. Covariance depends on the moments up to the order of the sum of the orders of the inputs, so in the present case it goes up to third order. $\endgroup$
    – Ben
    Aug 31, 2021 at 8:36
  • 1
    $\begingroup$ The sample skewness that you choose will certainly matter, since they are different and they have slightly different properties. However, all the major sample skewness statistics are consistent estimators of the true skewness, so regardless of which you choose you should get a consistent estimator. For comparisons of skewness estimators, see Joanes and Gill (1998). $\endgroup$
    – Ben
    Aug 31, 2021 at 8:38
  • $\begingroup$ Thanks very much for your comments, the intuition and the Joanes and Gill (1998) paper in particular!:) Perhaps one short final question about a possible generalization: Suppose there are two time series $\{x_t\}_{t=1,...,T}$ and $\{y_t\}_{t=1,...,T}$ with means $\mu_x$ and $\mu_y$ and variances $\sigma^2_x$ and $\sigma^2_y$. Can we also compute the $4\times4$ covariance matrix of $\hat{\mu_x}$, $\hat{\mu_y}$, $\hat{\sigma^2_x}$, and $\hat{\sigma^2_y}$? $\endgroup$
    – Alex
    Aug 31, 2021 at 9:41
  • $\begingroup$ That would be a fairly trivial extension if the two series are independent. $\endgroup$
    – Ben
    Aug 31, 2021 at 11:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.