There's more than one way to parameterize MA models. Your parameterization disagrees with the one used by the software.
Your model for the time series process $(Y_t)$ is
$$(Y_t - \mu) = (1 + \theta_1 L + \theta_2 L^2 + \theta_3 L^3 + \theta_4 L^4)\, \varepsilon_t$$
where $\mu=0$ is the mean, $L$ is the lag operator, the $\epsilon_t$ are white noise, and the parameters are $(\theta_i) = (1,0,0,-0.8).$
The four roots $\lambda_1, \ldots, \lambda_4$ of the associated polynomial $1 + \theta_1 z + \cdots + \theta_4 z^4 = 1 + z - 0.8 z^4$ are $-0.7486682, -0.2769154 \pm 1.09781i,$ and $1.3024989.$ Because they are not all outside the unit circle (the first one has a norm less than $1$), this process is not invertible.
There is an invertible process with exactly the same statistical properties, obtained by using the reciprocals of the roots inside the unit circle and rescaling the polynomial to have a constant term of $1.$ That replaces $-0.7486682$ by $-1.3357052$ and the associated polynomial becomes
$$\eqalign{
&(z - (-1.3357))(z - (-0.2769 + 1.0978i))((z - (-0.2769 - 1.0978i))(z - 1.3024) \\
&\propto 1 + 0.41296 z + 0.19707 z^2 - 0.26322 z^3 - 0.4484 z^4,
}$$
corresponding to the model
$$(Y_t - \mu) = (1 + 0.41296 L + 0.19707 L^2 - 0.26322 L^3 - 0.4484 L^4)\, \delta_t.$$
(The fundamental innovations $\delta_t$ are a constant multiple of the original white noise process $\varepsilon_t.$)
Because the software assumes the process is invertible, these are the coefficients it estimates. Indeed, when you don't constrain the estimates
> arima(x, order = c(0,0,r), method = "ML")
the output includes
Coefficients:
ma1 ma2 ma3 ma4 intercept
0.4238 0.1429 -0.2699 -0.4662 0.0173
s.e. 0.0402 0.0437 0.0449 0.0408 0.0493
and you can observe a close match of all estimates with the invertible coefficients (within $1.25$ times the standard errors in this case).
The solution to your quandary is to solve this exercise by proposing coefficients of an invertible model in the first place. You don't have to find polynomial roots to do that: simply pick $k$ numbers $\lambda_i$ of (complex) modulus no less than $1$ for the roots, construct the polynomial $$p(z) = \frac{\prod_{i=1}^k (z-\lambda_i) }{ \prod_{i=1}^k (-\lambda_i)} = 1 + p_1 z + p_2 z^2 + \cdots + p_k z^k,$$ and specify the vector $(p_1,p_2, \ldots, p_k)$ in the simulation.
Here's an example in R with $(\lambda_i)=(1,3/2,\pm 2i).$ Because the complex roots occur in conjugate pairs ($\pm 2i$), the model will be real:
lambda <- c(1, 3/2, 2i, -2i)
p <- 1; for (l in lambda) p <- convolve(p, c(1, -l), type="open")
(coeff <- Re(p[-1] / p[1]))
The coefficients are $(-5/3, 11/12, -5/12, 1/6) \approx (-1.667, 0.917, -0.417, 0.167).$ Now simulate and fit. (This time let's generate a much longer series so the estimates will be more reliable.)
set.seed(17)
x <- arima.sim(model = list(ma=coeff), n = 5000)
arima(x, order = c(0,0,length(coeff)), method = "ML")
Coefficients:
ma1 ma2 ma3 ma4 intercept
-1.6536 0.9069 -0.4110 0.1577 0
s.e. 0.0141 0.0268 0.0265 0.0141 0
sigma^2 estimated as 0.9754: log likelihood = -7037.65, aic = 14087.31
Not only are the coefficients well-estimated, so also is the variance of the innovations $\sigma^2$ (which was equal to $1$ by default).
Reference
Hamilton, James D. (1994), Time Series Analysis, Princeton University Press: Section 3.7, Invertibility.
