How does the the variance of a moving average process vary with its order? Let $\ MA(q)$ be a moving average process of order $\ q$, and $\hat{y} = \mu + a_t$ be a constant forecast model whose level is represented by $\mu $ and the random term at time $\ t$ by $\ a_t$. 
Now, suppose $\mu \equiv \hat{a_1}$ so one can write  $\hat{y} = \ a_1 + a_t$. 
If $\hat{a}_1$ is calculated by a $\ MA(q)$, what would the variance of $\hat{y}$ be? 
What if, then, $\hat{a}_1$ is calculated by a $\ MA(2q)$ or a $\ MA(q^2)$ - how would that affect the variance of $\hat{y}$?

My goal is understand how the variance of my forecasted values will vary as I vary$\ q$.

Bob Stine, from UPenn, made some actual calculations on page 9 of this document, but the way he calculates the variance for 3 steps ahead is something I'm yet to understand. More precisely, I don't get  why he keeps getting rid of his $ a_n$ terms. 
Thank you. 
..........................................................................................................................................................
Edit: I'll try to show some of my thoughts on this matter here: 
On page 86, Montgomery, Johnson,and Gardiner (1990) write that: 

$\large var(M_T) = \large var[\frac{1}{N}\sum_{k=0}^{N-1} x_{T-k}] = \large \frac{\sigma^2_\epsilon}{N}$

where $\ M_T$ is the $\ MA$ process, $\ N$ is the number of the most recent observations taken into account, and $\sigma^2_\epsilon$ is the variance of the noise (random component). 
However, if I write the $\ MA$ process as $\phi_p. z_t = \theta_q . a_t$, for $\phi_p$ equals 1, I have that my time series $z_t$ is equal to $\theta_q . a_t$, where $\theta_q$ is the lag operator of order $q$ and $a_t$ is my random component. 
Thus, I can write $z_t$ as $z_t = (1 + B + B^2 + ... + B^q).a_t$, and then $z_t = a_t + a_{t-1} + a_{t-2} + ... + a_{t-q+1}$. 
If then I take the variance of $z_t$, $var(z_t)$, I take the variance of $a_t + a_{t-1} + a_{t-2$} + ... + a_{t-q+1}$, which is the same as:
$var(z_t)$ = $\sum_{k=1}^{t-q+1} var(a_t)$. Now, if $var(a_t) = \sigma^2_\epsilon$, the variance of z will be equal to $N.\sigma^2_\epsilon$, where $N$ is defined as above. 
My result clearly doesn't match that of Montgomery et al. 
Assuming it's quite fair to say I made a mistake, where is it? 
 A: In order to keep it simple, I use a MA(q)-process where $\mu=0$ and $E[a_t] = 0$.

*

*Take the random shock model and express the MA(q) process as a linear filter, $y_t = \Psi(B) \cdot a_t = \sum_{j=0}^\infty \psi_j \cdot a_{t-j}$. Therefore,
\begin{align}
y_{t+\ell} &=  \sum_{j=0}^\infty \psi_j \cdot a_{t + \ell -j} \\
&= \sum_{j=0}^{\ell-1} ... + \sum_{j=\ell}^{\infty} ...
\end{align}
Looking at this eq. from the time $t$ (in the forward direction), the first sum contains only future terms, while the second contains only past terms.

*Suppose that the "best" forcast (smallest RMS) at time $t$ and lead time $\ell$ is given by
$$\hat{y}_t(\ell) = \sum_{j=0}^\infty \psi^\star_{j+\ell} a_{t-j}$$
where the $\psi^\star_{j+\ell}$ need to be determined. Note, that the forecast contains only past terms. Hence, using the eq. from (1) we know that
$$E[y_{t+\ell} - \hat{y}_t(\ell)]^2 = 
\sum_{j=0}^{\ell-1} \psi^2_j \cdot \sigma^2_a 
+ \sum_{j=\ell}^{\infty} (\psi_{\ell+j} - \psi^\star_{\ell + j} )^2 \sigma^2_a$$
where the first term contains the future terms and the second the past terms. Clearly this is minimized by setting $\psi^\star_{\ell + j} = \psi_{\ell+j}$.

*Put this result back to eq. (1) we identify the to sums,
\begin{align}
y_{t+\ell} &=  \sum_{j=0}^{\ell-1} ... + \sum_{j=\ell}^{\infty} ... \\
&= e_t(\ell) + \hat{y}_t(\ell)
\end{align}
as the error of the forecast, $e_t(\ell)$, and the forcast, $\hat{y}_t(\ell)$, at time $t$ and lead time $\ell$, respectively.

*Now that we have an expression for the error of the forecast, we use the index of $a_t$ to obtain its variance,
$$var[e_{t}(\ell)] = \sum_{j=0}^{\ell-1} \psi^2_j \sigma_a^2$$
Finally, note that this calc is valid for an ARMA(p, q)-process. For a MA(q) we simply use $\psi_j =0 $ for $j>q$.
