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.]

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.
