Variance of sum of AR(2) processes We have
$$y_t = a\mu_t + b\mu_{t-1} + c\mu_{t-2} + d\mu_{t-3}$$
and $\mu$ itself is an AR process, let's assume here an AR(2)
$$\mu_t = \phi_1 \mu_{t-1} + \phi_2 \mu_{t-2} + \epsilon_t$$
where $\epsilon_t$ is a white noise process with E($\epsilon_t$)=0 and $Var(\epsilon_t)=\sigma^2$.
As there will be covariance terms present, how can one solve (or compute numerically) the variance of $y_t$? 
Note that this question is close to this one (which is still unanswered); the the goal here is however to find an approach that excludes any covariance terms between different lags of the $\mu$s, as the sum in the $y_t$ equations could get long. 
My approach so far 


*

*Rewrite the AR process in companion form:
$$
\begin{bmatrix}
\mu_{t} \\
\mu_{t-1} 
\end{bmatrix}=
\begin{bmatrix}
\phi_1 & \phi_2 \\
1 & 0 
\end{bmatrix}
\begin{bmatrix}
\mu_{t-1} \\
\mu_{t-2} 
\end{bmatrix}+
\begin{bmatrix}
\epsilon_{t} \\
0
\end{bmatrix}
$$
$$\pmb \mu_t = \pmb \Phi \pmb \mu_{t-1} + \pmb \epsilon_t$$


*Normally, I would now compute the variance of the $\mu$s by a solver for Lyapunov equations (in either R or Matlab). (Rewriting the AR as $\Sigma = \Phi \Sigma \Phi' + Q$, and then solving for $\Sigma$). Yet here, we have to deal with the covariances between the lags. 


Ideally, I would now sub out the lags of the $\mu$s in the last equation, as in an AR(1). But can I do this? 
$$\pmb \mu_t = \sum_{i=0} \pmb \Phi^i \pmb \epsilon_{t-i}$$


*Once the $\mu$s are expressed in terms of $\epsilon$s only, one could calculate the variances of each $\mu$ term in the $y_t$ equation separately and ignore any covariances


Please note:
The ultimate goal is to estimate $Var(y_t)$; parameters $a,b,c,d,\Phi$ and $\sigma$ are known.
 A: Writing the model using backshift operator polynomials, you have
$$
y_t = \theta_3(B)\mu_t \tag{1}
$$
and
$$
\phi_2(B)\mu_t = \epsilon_t. \tag{2}
$$
Applying $\phi_2(B)$ to both sides of (1) and using (2) yields
$$
\phi_2(B)y_t = \theta_3(B)\phi_2(B)\mu_t = \theta_3(B)\epsilon_t
$$
which shows that $y_t$ is ARMA(2,3) with autocovariance function that can be computed using standard methods, implemented in e.g. R-function ltsa:::tacvfARMA.  
Note that if $a\neq 1$, you'll need to do some rescaling to make your model conform to the convention that the first coefficients of the the MA and AR-polynomial are 1.
A: Evidently $\operatorname{Var}(y_t)$ will involve covariances of the $\mu$ series at lags $0,$ $1,$ $2,$ and $3.$  A comment to the question you reference directs us to a thread that explains how to find those covariances.  You should obtain the solution
$$\eqalign{
&\operatorname{Var}(\mu_t,\mu_t) &= \gamma_0 = \left(\frac{1-\phi_2}{1+\phi_2}\right)\frac{\sigma^2}{(1-\phi_2)^2-\phi_1^2} \\
&\operatorname{Var}(\mu_t,\mu_{t-1}) &=\gamma_1 = \frac{\phi_1}{1-\phi_2}\,\gamma_0 \\
&\operatorname{Var}(\mu_t,\mu_{t-2}) &=\gamma_2 = \phi_1 \gamma_1 + \phi_2 \gamma_0 \\
&\operatorname{Var}(\mu_t,\mu_{t-3}) &=\gamma_3 = \phi_1 \gamma_2 + \phi_2 \gamma_1.
}$$
Because covariance is bilinear,
$$\eqalign{
\operatorname{Var}(y_t) &= \operatorname{Var}((a,b,c,d)(\mu_t,\mu_{t-1},\mu_{t-2},\mu_{t-3})^\prime) \\
&= \pmatrix{a&b&c&d}\Gamma \pmatrix{a\\b\\c\\d}
}$$
where
$$\Gamma = (\gamma_{|i-j|}) = \pmatrix{\gamma_0 & \gamma_1 & \gamma_2 & \gamma_3 \\
\gamma_1 & \gamma_0 &\gamma_1 &\gamma_2 \\
\gamma_2 & \gamma_1 & \gamma_0 & \gamma_1 \\
\gamma_3 & \gamma_2 & \gamma_1 & \gamma_0}.
$$
