Arima simulation followed by modeling in R produces bad estimates

Simulate a moving average using arima.sim in R. Then estimate the coefficients using arima. Fix all coefficients to their true values, except one:

set.seed(222)
x <- arima.sim(model = list(ma=c(1, 0, 0, -0.8)), n = 500, rand.gen = rnorm)
arima(x,order = c(0,0,4), fixed = c(1, 0, 0, NA, 0), method = "ML")


Results:

Call:
arima(x = x, order = c(0, 0, 4), fixed = c(1, 0, 0, NA, 0), method = "ML")

Coefficients:
ma1  ma2  ma3      ma4  intercept
1    0    0  -3.6772          0
s.e.    0    0    0   0.2970          0

sigma^2 estimated as 0.159:  log likelihood = -901.13,  aic = 1806.26


I expect the estimate of ma4 to be close to $$-0.8$$, but get $$-3.68$$. Why is the parameter estimate so far from the true parameter?

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.