Covariance of the means of $x_t$ and $y_t=(x_t-\bar{x})^2$ 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}$ 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.
