Autocorrelation of a MA(1) process with correlated noise sources I have the following problem. I want to find the expression for the autocorrelation of the following random variable or MA(1) process.
$X_n = \epsilon_n - \epsilon_{n-1}$
Where $\epsilon$ is a stationary correlated noise with finite mean $\mu_\epsilon$ and standard deviation $\sigma_\epsilon$.
So the mean and variance of $X$ is $\mu_X = 0$ and $\sigma_X^2 = 2\sigma^2_\epsilon$ (this was the main mistake in the derivation)
I start from this expression for the autocorrelation.
$ \displaystyle  R^k_X = \frac{E[(X_n-\mu_X) (X_{n+k}-\mu_X)]}{\sigma_X^2} $
For simplicity I worked with the lag-one, so this is:
$ \displaystyle  R^1_X = \frac{E[(\epsilon_n-\epsilon_{n-1})(\epsilon_{n+1}-\epsilon_{n})]}{2\sigma^2_\epsilon}$
Expanding and taking the expectancy it gives:
$ \displaystyle  R^1_X = \frac{E[\epsilon_n\epsilon_{n+1}] +
 E[\epsilon_n\epsilon_{n-1}]-E[\epsilon_{n+1}\epsilon_{n-1}]-E[\epsilon^2_n]}{2\sigma^2_\epsilon}$
The first two terms are equivalent (correct me if I am wrong), and they are related to the lag one autocovariance of the noise, but they don't have their means subtracted, and the third term is related to the lag two. Using $E[\epsilon_n \epsilon_{n+k}] = Cov[\epsilon_n, \epsilon_{n+k}] + \mu^2_\epsilon$ and adding and subtracting the squared means allow rewriting to:
$ \displaystyle  R^1_X = R^1_\epsilon -\frac{R^2_\epsilon}{2} + \frac{\mu_\epsilon^2 - E[\epsilon^2_n]}{2\sigma^2_\epsilon}$
Then for the last term, intead of covariance we have variance $E[\epsilon_n^2] = \sigma^2_\epsilon + \mu^2_\epsilon$ and this gives:
$ \displaystyle  R^1_X = R^1_\epsilon - \frac{R^2_\epsilon}{2} - 0.5$
The problem arises when checking this expression with simulated processes. For the correlated $\epsilon$ noise I used another MA(1) process with uncorrelated normal noise source $\xi$, thus having something like $\epsilon_n = \xi_n+\alpha \xi_{n-1}$, so the lag-two term is theoretically zero. 
The following picture shows the deviations from the value of -0.5 for different $\alpha$, it is only valid when $\alpha$ is 0 or 1. I hope there is a problem in the deviation of the formula.

This is the python/numpy code for the simulation, R1 measures the lag-one crosscorrelation:
def R1(x,y):
    return mean((x[1:]-x.mean())*(y[:-1]-y.mean()))/std(x)/std(y)

xi = randn(100000)

eps = xi[1:] + alpha*xi[:-1]
X = eps[1:] - eps[:-1]

out = R1(X,X) - R1(eps,eps)

 A: Theory
Adding a constant (say $\lambda$) to $\epsilon_n$ makes no difference, because $$(\epsilon_n + \lambda) - (\epsilon_{n-1} + \lambda) = \epsilon_n - \epsilon_{n-1} = x_n.$$  By choosing $\lambda=-\mu_\epsilon$ you may assume $E[\epsilon_n]=0$.
This makes covariances easy to compute in terms of the lagged correlations of $(\epsilon_n)$ given by $E[\epsilon_{n+i}\epsilon_n]=\rho_i \sigma_\epsilon^2$ for any $i$:
$$\eqalign{\operatorname{Cov}(x_{n+j}, x_{n}) &= E[x_{n+j}x_{n}]= E[(\epsilon_{n+j}-\epsilon_{n+j-1})(\epsilon_{n}-\epsilon_{n-1})]\\&=E[\epsilon_{n+j}\epsilon_n] - E[\epsilon_{n+j}\epsilon_{n-1}] - E[\epsilon_{n+j-1}\epsilon_n] + E[\epsilon_{n+j-1}\epsilon_{n-1}] \\&= \sigma_\epsilon^2(\rho_j - \rho_{j+1} - \rho_{j-1} + \rho_j) = -\sigma_\epsilon^2(\rho_{j-1} - 2\rho_{j} + \rho_{j+1}).}$$
A particular case of this is the variance, whose expression can be simplified from the facts $\rho_{-i}=\rho_i$ and $\rho_0=1$:
$$\operatorname{Var}(x_n) = \operatorname{Cov}(x_n,x_n) = -\sigma_\epsilon^2(\rho_{-1} - 2\rho_0 + \rho_1)=2\sigma_\epsilon^2(1-\rho_{1}).$$
Dividing the covariances by the variance  yields the lag-$j$ correlation of $(x_n)$, written $\rho_j(x)$, in terms of the correlations of $(\epsilon_n)$,
$$\rho_j(x) = -\frac{\rho_{j-1} - 2\rho_{j} + \rho_{j+1}}{2(1-\rho_1)}.\tag{1}$$
Similar calculations from the definition $\epsilon_n = \xi_n + \alpha \xi_{n-1}$ give
$$\rho_1(\epsilon) = \frac{\alpha}{1+\alpha^2}\tag{2}$$
and $\rho_j(\epsilon) =0$ for $j\gt 1$.  Plugging $(2)$ into $1$ and simplifying gives
$$\rho_1(x) = -\frac{1}{2}\frac{1-2\alpha+\alpha^2}{1-\alpha+\alpha^2}.$$
Note that this implies $-2/3 \le \rho_1(x) \le 0$, a range that is fully covered as $\alpha$ ranges from $-1$ through $1$.
Simulation
$50$ separate simulations of $(\xi_n)$ were used to produce $21$ simulated series $(x_n)$, $n=1,\ldots, 10,000$, for an equally spaced range of values of $\rho_1$ from $-2/3$ through $0$. (The R code is at the end of this message.)  The plots of the observed lag-1 autocorrelation coefficients against the values of $\rho_1$ (shown below in red) form a tight envelope around the line of equality (shown below in gray), visibly demonstrating the correctness of these calculations.

We must conclude that the code used to create the plot in the question is erroneous.  
n <- 1e4                            # Length of (x) series
r <- seq(-2/3, 0, length.out=21)    # Array of rho_1(x) values to study
alpha <- ifelse(r==-1/2, 0, (1 + r - sqrt(-2*r - 3*r^2)) / (1 + 2*r))
#
# Simulate (x) many times.
#
X <- lapply(1:50, function(i) {
  xi <- rnorm(n+2)
  sapply(alpha, function(a) {
    epsilon <- xi[-1] + xi[-(n+2)] * a
    x <- epsilon[-(n+1)] - epsilon[-1]
    acf(ts(x), lag.max=1, plot=FALSE)$acf[2,1,1]
    # acf(ts(epsilon), lag.max=2, plot=FALSE)$acf[2,1,1]
  })
})
#
# Prepare to plot the simulated results.
#
rho.1 <- alpha / (1 + alpha^2)    # Lag-1 acf of epsilon
rho.2 <- 0                        # Lag-2 acf of epsilon
rho <- (2*rho.1 - (1 + rho.2)) / (2*(1-rho.1))    # Lag-1 acf of x
#
# Draw a plot window and reference line.
#
plot(range(rho), range(X), type="n", asp=1,
     xlab=expression(rho[1]), ylab="Observed",
     main="Lag-1 Correlations of X")
abline(c(0,1), col="gray", lwd=2)
#
# Plot the results.
# (Plotting them as points is commented out.)
# invisible(lapply(X, function(y) points(rho, y, pch=16, col="#ff000010")))
invisible(lapply(X, function(y) lines(rho, y, pch=16, col="#ff000020")))

