Estimating the error in the average of correlated values tl;dr I can only generate samples of a random variable $X$ using MCMC.  How can I find the error in the estimate of the expected value of $X$ based on this MCMC data?

The problem
I have a "black box" which gives me a sequence of "random" numbers (MCMC sampler).  I have recorded $n$ numbers from this black box, which is all the data I can work with (I can't sample additional points).  Now I need to estimate the average of the black box based on this data, and get some idea about how precise the estimate is.
In other words, if I simply calculate the average of these $n$ datapoints, how far is it going to be from the true average of the black box?
The difficulty: Consecutive numbers coming form the black box are strongly correlated, however the $i^\text{th}$ and $(i+k)^\text{th}$ numbers are effectively uncorrelated if $k$ is large.
Are there any methods to deal with this type of situation, and find the error in the estimate?  Let's assume that $n$ is large enough that I have some uncorreltaed data points.

What I'm doing now
In the following I'll describe what I am doing right now and why I am not satisfied with this method.  This method works but it is not ideal.
Let's look at the case when the numbers $x$ are uncorrelated.  Let $\langle x \rangle_n$ denote the average of $n$ consecutive numbers, which is an estimate of the expected value of the black box.  Then $\operatorname{Var}(\langle x \rangle_n) = \operatorname{Var}(x)/n$.
$\sqrt{\operatorname{Var}(\langle x \rangle_n)}$ can be used as the error of the estimate $\langle x \rangle_n)$ and it's easy to compute based on $\operatorname{Var}(x)$.
However my actual data is correlated.  So I partitioned my dataset of $n$ elements into runs of $k$ elements and calculated the product $k \times \operatorname{Var}(\langle x \rangle_k)$ for various $k$.  When $k$ is large enough, and the correlations are not important any more, this product will settle at a constant value, which can be used to extrapolate for $\operatorname{Var}(\langle x \rangle_n)$.  This is a plot of $k \operatorname{Var}(\langle x \rangle_k)$ for various values of $k$ for four different datasets (four colours) of size $n = 2^{27}$:

You can see that as $k$ becomes large enough, the product $k \operatorname{Var}(\langle x \rangle_k)$ settles on a constant.  But as $k$ becomes too large, I can only parition my length-$n$ dataset into a few parts of length $k$, so the estimate of the variance becomes inaccurate.  Hence the noise on the right of the curve.
So I need to find the flat part of the curve and discard the noisy part, then estimate the variance of the average of the full dataset ($\langle x \rangle_n$) based on it.
Why my method isn't ideal
The method I describe above is complicated and involves several steps, and several choices in each step.  For how many (and what magnitude of) $k$ values should I compute the variances?  How do I determine where is the flat part of the curve?  How do I decide which part is too noisy to keep?  All of these are separate problems on their own with many valid heuristic solutions, which will all influence the result.
Also, this method isn't very efficient computationally.
I am hoping for a simpler and cleaner robust method where the details won't have such a large effect on the end result.
 A: This is an active area of research. This first question is whether a central limit theorem (CLT) even exists which depends on the mixing properties of your MCMC, e.g. geometric convergence. Typically, this is a nontrivial question. 
Provided a CLT exists, the second question is how to obtain a consistent estimator of the variance in the CLT. Partitioning the chain as you have done is called batching and is the approach used in the R package mcmcse. I'd suggest you take a look at that package and the references therein. As an example of using the package you can do 
library(mcmcse)
n = 1000
y = rep(0,n)
for (i in 2:n) y[i] = rnorm(1,0.9*y[i-1])
mcse(y)

which will return an estimate of the expectation of y, which is known to be zero in this case, as well as an estimate of the Monte Carlo standard error in the estimator which takes into account the correlation in the chain.
This package allows the user to specify the batch size, but defaults to the square root of the total number of samples. Perhaps the references in the package will give an indication of why this is the default or what other batch sizes will work.
A: The problem of finding error estimates of statistics in (autocorrelated) time series is usually approached via block bootstrapping. It is the same in spirit as your approach. See Section 5 of this document for a very short summary [1]. There is also some parallel work in the physics community, where ideas from renormalisation are used; see e.g. [2] for a clear exposition.
References:
[1] Kreiss, J. P., & Lahiri, S. N. (2012). Bootstrap methods for time series. In Handbook of statistics (Vol. 30, pp. 3-26). Elsevier.
[2] Flyvbjerg, H., & Petersen, H. G. (1989). Error estimates on averages of correlated data. The Journal of Chemical Physics, 91(1), 461-466.
