arima.sim (again) what does the mean parameter do in an ma(1) sim

Question

When I run the following sim:

set.seed(999)
inter.vec <- vector()
mean.vec <- vector()
for (i in 1:999){
    yt <- arima.sim(n=9999, list(order=c(0,0,0)), mean = 0.6)
    mean.vec[i] <- mean(yt)

    fit <- Arima(yt, order = c(0,0,0), include.mean = T)
    inter.vec[i] <- coef(fit)[1]
}
lattice::densityplot(mean.vec)
lattice::densityplot(inter.vec)

I get exactly what I expect under the specification $Y_t = \mu + e_t$. The means of the generated series are centred on 0.6 and the estimate for the intercept from the Arima model is centred on 0.6.

Then I run the following sim, an MA(1) model with non-zero mean, initially under the impression I am generating data from the specification $Y_t = 0.6 -0.3 e_{t-1} + e_t$:

set.seed(616)
ma.vec <- vector()
inter.vec <- vector()
mean.vec <- vector()
for (i in 1:999){
   yt <- arima.sim(n=9999, list(order=c(0,0,1), ma=-.3), mean = 0.6)
   mean.vec[i] <- mean(yt)

   fit <- Arima(yt, order = c(0,0,1), include.mean = T)
   ma.vec[i] <- coef(fit)[1]
   inter.vec[i] <- coef(fit)[2]
}

But the model that I thought I was generating data from is not recovered as can be seen from the following plots. Clearly, the MA parameter is recovered well, but the intercept/mean is not. So what model is the sim generating data from?

Answer 1

For an MA(1) models (as per my question), the implementation of arima.sim initially constructs a sample of random values (innov in the function call) using the mean value that I passed in ($0.6$). If the user does not say otherwise, a normal distribution is assumed giving a series of random values distributed as $e_i \sim N(0.6, 1)$. As the sample size increases the mean of the realised series will tend to 0.6, i.e. $\bar{e} = \frac{\sum_i^n e_i}{n} \approx 0.6$.

Next, arima.sim makes use of the stats::filter function, which in my case is passed a filter parameter equal to the vector c(1, -0.3). In response, the filter then sums the current value and the product of -0.3 and the previous value and this is the requested MA(1) series.

However, the new mean will be as follows:

\begin{align} \bar{e}_{new} &= \frac{(-0.3e_1 + e_2)+(-0.3e_2 + e_3)+...+(-0.3e_{n-1} + e_n)}{n} \\ &= \frac{-0.3e_1 + 0.7e_2 + 0.7e_3+...+0.7e_{n-1} + e_n}{n} \end{align}

For simplicity apply $E[e_1] = E[e_n] = 0$ gives:

\begin{align} \bar{e}_{new} &= \frac{0.7e_2 + 0.7e_3+...+0.7e_{n-1}}{n}\\ & \approx 0.7 \bar{e} \approx 0.7 \times 0.6 = 0.42 \end{align}

So, it looks to me like the underlying data generating process is (approximately) \begin{align} Y_t=0.42 - 0.3 e_{t-1}+e_t \end{align} and not $Y_t=0.6 - 0.3 e_{t-1}+e_t$ as one may be (possibly) forgiven for thinking.