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.]
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:
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:
and similarly the OP's sample variance estimator for $Y$ in terms of bivariate power sums is:
The list of our 4 estimators is then:
Then the variance-covariance matrix of the 4 estimators is:
where:
$$\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.
