The way an MA(q) model works I am trying to understand the way MA(q) models work.
For this purpose I have created a simple data set with only
three values. I then adapted a MA(1) model to it. The results
are shown below:
x<-c(2,5,3)
m<-arima(x,order=c(0,0,1))

Series: x 
ARIMA(0,0,1) with non-zero mean 

Coefficients:
          ma1  intercept
      -1.0000     3.5000
s.e.   0.8165     0.3163

sigma^2 estimated as 0.5:  log likelihood=-3.91
AIC=13.82   AICc=-10.18   BIC=11.12

While the MA(1) model looks like this: 
$$X_t = c +a_t - \theta*a_{t-1}$$
and $a_t$ is White Noise.
What I cant figure out is how to get the fitted values:
library(forecast)
fitted(m)
Time Series:
Start = 1 
End = 3 
Frequency = 1 
[1] 3.060660 4.387627 3.000000

I tried different ways, but I cant find out how the fitted values (3.060660, 4.387627 and 3.000000) are calculated.
I would be very thankful for an answer!
 A: I found it out, but I can't tell you the exact reason for it:
The problem is, that the initialization is unkown/strange, if you do an example with more values, you will see, that a simple MA(1) forecasting according to the following rule will work (notation in R of the MA is slightly different to yours, the sign of the theta is different):
\begin{align}\hat{X}_{T|T-1}=E(c+a_T+\theta∗a_{T−1})=c+\theta*a_{T-1}\end{align}
You can calculate these values manually, consider the following example:
z<-c(2,5,3,4,3,4,5,4.3,4.3,4.5,4.3,4.5,3.4,5.3,4.2,3.4,2.3,2.3,4.5,3.4,5,5.4,5.4,3.4,5.43,5.64,5.6,3.4,5.3,5,6.3,4.5)
m<-arima(z,order=c(0,0,1))

This gives
Series: z 
ARIMA(0,0,1) with non-zero mean 

Coefficients:
         ma1  intercept
      0.1162     4.2748
s.e.  0.1500     0.2076

sigma^2 estimated as 1.11:  log likelihood=-47.09
AIC=100.18   AICc=101.03   BIC=104.57

The values of the output can be used via
m$coef[1] and m$coef[1]
So you now compare the values with the following code:
m$coef[2]+ m$coef[1]*(2-4.259580)
fitted(m)[2]
m$coef[2]+ m$coef[1]*(5-4.013978)
fitted(m)[3]
m$coef[2]+ m$coef[1]*(3-4.389402)
fitted(m)[4]
m$coef[2]+ m$coef[1]*(4-4.113285)
fitted(m)[5]
m$coef[2]+ m$coef[1]*(3-4.261626)
fitted(m)[6]

You will notice, that at the beginning, there is a small difference, at the end, the values are the same! So the simple forecasting rule of a MA(1) does hold, but R seems to do some specific initialization. I know that STATA uses a certain Kalman filter setting, maybe R is doing the same. I hope this helps.
If it did help you, you can accept my answer by clicking on the hook to the left of my post.
A: In terms of where those particular numbers come from, this seems to do the trick:
> x - m$residuals
Time Series:
Start = 1 
End = 3 
Frequency = 1 
[1] 3.060660 4.387627 3.000000

Or
> x + (m$residuals * ma1)
Time Series:
Start = 1 
End = 3 
Frequency = 1 
[1] 3.060660 4.387627 3.000000

I don't know in what sense these are meant to be the 'fitted' values (maybe someone else can help?). For an ARIMA model our explanatory data are also our fitted values. 
Our estimator for $Y_k$ given the history up to time $T$, $H_T$, where $T>k$, is $\hat y_k = \mathbb E [Y_k | H_T ] = y_k$, i.e. we already have the answer in the process history.  So we have to ask, our fitted values anchored at what time? For a length $T$ series I guess there would be $T$ possible fitted values for each point. 
If anyone asks, based on having the data $y_{1:T}$ what's the best estimate of $y_k$ - seems there's only one pragmatic answer, $y_k$. The residuals are interesting in the sense that if they're massive relative to the data, that's saying something about the magnitude of the process variance.
