You could obtain the fitted values as one-step forecasts using the innovations algorithm. See for example proposition 5.5.2 in Brockwell and Davis; downloable from the internet I found these slides.
It is much easier to obtain the fitted values as the difference between the observed values and the residuals. In this case, your question boils down to
obtaining the residuals.
Let's take this series generated as an MA(1) process:
set.seed(123)
x <- arima.sim(n = 150, model = list(ma = 0.4))
fit <- arima(x, order = c(0,0,1), include.mean = TRUE)
resid <- residuals(fit)
The residuals, $\hat{e}_t$, can be obtained as a recursive filter:
$$
\hat{e}_t = x_t - \hat{\mu} - \hat{\theta} \hat{e}_{t-1}
$$
For example, we can obtain the residual at time point $140$ as the observed value at $t=140$ minus the estimated mean minus $\hat{\theta}$ times the previous residual, $t=139$):
macoef <- coef(fit)[1]
mu <- coef(fit)[2]
as.vector(x[140] - mu - macoef * resid[139])
# [1] 0.7742192
# equal to
residuals(fit)[140]
# [1] 0.7742192
The function filter can be used to do these calculations:
resid.v2 <- filter(x = x - mu, filter = -macoef, method = "recursive")
You can see that the result are very close to the residuals returned by residuals. The difference in the first residuals is most likely due to some initialization that I may have omitted.
head(cbind(resid, resid.v2))
# resid resid.v2
# [1,] -0.39447063 -0.429102263
# [2,] 1.62425953 1.675613161
# [3,] 0.03344943 0.001881606
# [4,] 0.16839438 0.181951010
# [5,] 1.71983927 1.714147387
# [6,] 0.43595315 0.438334539
tail(cbind(resid, resid.v2))
# resid resid.v2
# [145,] -0.4803322 -0.4803322
# [146,] -1.4432094 -1.4432094
# [147,] 0.7463573 0.7463573
# [148,] 2.0810053 2.0810053
# [149,] -1.3126564 -1.3126564
# [150,] 0.8601761 0.8601761
The fitted values are just the observed values minus the residuals:
require(forecast)
head(cbind(x - resid.v2, fitted(fit)))
# [1,] -0.02526549 -0.05989712
# [2,] -0.20897584 -0.15762221
# [3,] 0.69211011 0.66054229
# [4,] -0.02445992 -0.01090329
# [5,] 0.05263269 0.04694081
# [6,] 0.70860766 0.71098905
tail(cbind(x - resid.v2, fitted(fit)))
# [145,] -0.6911888 -0.6911888
# [146,] -0.2309088 -0.2309088
# [147,] -0.6431427 -0.6431427
# [148,] 0.2942704 0.2942704
# [149,] 0.8656695 0.8656695
# [150,] -0.5872494 -0.5872494
In practice you should use the functions residuals and fitted but for pedagogical purpose you can try the recursive equation used above. You can start by doing some examples by hand as shown above. I recommend you to read also the documentation of function filter and compare some of your calculations with it. Once you understand the operations involved in the computation of the residuals and fitted values you will be able to make a knowledgeable use of the more practical functions residuals and fitted.
You may find some other information related to your question in this post.
