Estimating unconditional variance in time series Consider a time series process with a well-defined, finite unconditional variance. Given a realization of the process (a time series) and a model for it, there are at least two ways of estimating the unconditional variance:


*

*Use the model to derive an analytical expression of the unconditional variance as a function of error variance and model parameters, and substitute sample estimates for the theoretical quantities.*

*Use sample variance as if the data were $i.i.d.$**


Question: When should 1. be preferred over 2. and vice versa?
A partly overlapping question is this one.

I simulated some AR(1) processes and discovered that 1. and 2. give very similar results. The kernel densities of the estimated variances over 1000 repetitions are nearly identical, and this holds across different sample sizes and different values of the autoregressive coefficient.
I then simulated some GARCH(1,1) processes and found that 2. is a much more robust alternative. While the median estimates due to 1. and 2. over 1000 repetitions were almost the same, 1. had many more extreme values that were apparently due to occasionally poor estimation of the GARCH model parameters.
So far it looks as if the simpler alternative 2. fares quite well compared with a more sophisticated alternative 1...

*Example: for an AR(1) process, the unconditional variance is $\frac{\sigma^2}{1-\varphi_1^2}$ where $\sigma^2$ is the error variance and $\varphi_1$ is the autoregressive coefficient. To make the estimator operational, $\sigma^2$ and $\varphi_1$ have to be substituted by their sample counterparts due to, say, full maximum likelihood estimator for the AR(1) model. 
** That is, $\widehat{\text{Var}}(X)=\frac{1}{n-1}\sum_{t=1}^T (x_t-\bar{x})^2$.
 A: The more population parameters one has to estimate, the worse. Moreover, the model may be misspecified.
On the other hand, the sample variance does not depend on the assumption of the elements of the process being independent, for its good asymptotic properties. We assume that these elements are homoskedastic, and this makes the sample variance a consistent estimator, irrespective of whether the elements are autocorrelated. 
More formally, assume $\{X\}$ is an ergodic process for its mean, so $E(X)$ is constant and finite, and the sample mean of a realization of $\{X\}$ is a consistent estimator for $E(X)$ (by Kinchine's Theorem). By the properties of ergodicitiy, $\{Z\}\equiv \{X^2\}$ is also an ergodic process for its mean, and $W\equiv Z+b$ is also an ergodic process for its mean, for $b$ a constant. Set $b\equiv -(E[X])^2$. Then
$$E(W) = E(Z) -(E[X])^2 = E(X^2) -(E[X])^2=\text{Var}(X)$$
So $E(W)$ is consistently estimated by its sample counterpart, which is the sample variance of $\{X\}$.
A: One solution would be to use both estimators simultaneously in a GMM framework.  This would have the benefit of minimizing your standard error beyond that of either estimator used alone.  Furthermore, because the GMM model would be over-identified, you could test whether or not the two estimators are significantly different.  If they are, it would suggest that you choose between the estimators* or that your model used in estimating (1) does not adequately fit the data.  
At the very least, this will give you some insight into the stability of the estimators on a case by case basis.  
*Your simulations would suggest that you choose estimator (2).  This is ultimately a subjective decision, but given the guaranteed consistency of (2) and your simulations it would seem to be the reasonable choice.
GMM Estimator
Let $\theta_0$ be the true unconditional variance you wish to estimate.  The two estimators you proposed in respective order are $\hat \theta_1$ and $\hat \theta_2$.  You can obtain a third, more efficient, estimator $\theta^\dagger$ by minimizing the following quadratic:
$$
\theta^\dagger=\min_{\theta}\bigg\{\begin{pmatrix}
\hat \theta_1-\theta\\
\hat \theta_2-\theta 
\end{pmatrix}^T 
\hat W
\begin{pmatrix}
\hat \theta_1-\theta\\
\hat \theta_2-\theta 
\end{pmatrix}\bigg\}
$$
Where $\hat W$ is the estimated inverse covariance matrix of $(\hat \theta_1,\hat \theta_2)$.  Due to the non iid nature of the data, the variance of $\hat \theta_2$ (and the covariance) will have to be estimated with special care, whether that be through a block bootstrap technique, HAC standard errors, Newey-West, or whatever else you deem appropriate. 
Over-Identification and the J-Test
Because the GMM estimate is over-identified you could perform a "J-test" with the following hypotheses:
$$
H_0: \begin{pmatrix}
\theta_1-\theta_0\\
\theta_2-\theta_0 
\end{pmatrix} = \mathbf{0}\;\;\;\;\;\;
H_1: \begin{pmatrix}
\theta_1-\theta_0\\
\theta_2-\theta_0 
\end{pmatrix} \neq \mathbf{0}
$$
The test statistic would be $$J=\begin{pmatrix}
\hat \theta_1-\theta^\dagger\\
\hat \theta_2-\theta^\dagger
\end{pmatrix}^T 
\hat W
\begin{pmatrix}
\hat \theta_1-\theta^\dagger\\
\hat \theta_2-\theta^\dagger 
\end{pmatrix} \stackrel{d}{\rightarrow} \chi^2_1$$
which, as noted, has an asymptotic chi-squared distribution with 1 degree of freedom under the null.
Under the null hypothesis: $\theta^\dagger$ would be the most efficient estimator for the unconditional variance.  
Under the alternative: there is a discrepancy between $\hat \theta_1$ and $\hat \theta_2$.  Either one estimator is unstable or the model chosen to estimate $\hat \theta_1$ does not fit the data well.  Given the findings of your simulations, this would suggest that you choose $\hat \theta_2$ to estimate the variance or choose a different model to obtain $\hat \theta_1$.
Related Tests
Another related test that you may consider using is the Hausman-Wu test.  This test will only work if the variance of $\hat \theta_2$ is greater than that of $\hat \theta_1$.  In this setting it can be used to test the consistency of $\hat \theta_1$.  Rather than explain the test in full detail, I refer you to the Wikipedia page which explains the test statistic and asymptotic distribution. 
I know this doesn't directly answer your question, but at least it gives you a way to check things and perhaps answer the question yourself on a case by case basis.  In some cases, it may even help you form an improved, more efficient estimator.
