How to compute the statistics of the DFT of correlated noise?
After many weeks I give the answer to my own question.

There is a limit in which we can solve this problem in a reasonably simple way. Suppose we sum enough points in our DFT that the central limit theorem guarantees that the distribution of the sum's real and imaginary parts are Gaussian distributed. Then we only need to compute the variance. If we specialize the case of the variances of the real part of the DFT we can write

$$
\langle (\textrm{Re} X_k) (\textrm{Re} X_l) \rangle = \frac{1}{N^2}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}\langle x_n x_m \rangle \cos(2\pi n k / N)\cos(2 \pi m l / N)
$$
The thing to note is that $\langle x_n x_m \rangle$ is just the correlation function of the time domain samples, which we can denote $\rho(n-m)$. Putting this in, we get

$$
\langle (\textrm{Re} X_k) (\textrm{Re} X_l) \rangle = \frac{1}{N^2}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1} \rho(n-m) \cos(2\pi n k / N)\cos(2 \pi m l / N)
$$

This form in terms of the correlation function is useful because the correlation function is frequently known; By the Wiener-Khinchin theorem the correlation function is the Fourier transform of the spectral density.

This sum can be computed numerically, or even analytically for some particular forms of $\rho$. To compute covariance with imaginary parts of $X$ just put $\sin$ instead of $\cos$.

I found the idea for this in J. Schoukens and J. Renneboog, IEEE Transactions on Instrumentation and Measurement, Vol. IM-35, No. 3, September (1986).