Autocovariance of an ARMA(2,1) process - derivation of analytical model for $\gamma( k)$ I need to derive analytic expressions for the autocovariance function $\gamma\left(k\right)$ of an ARMA(2,1) process denoted by:
$y_t=\phi_1y_{t-1}+\phi_2y_{t-2}+\theta_1\epsilon_{t-1}+\epsilon_t$
So, I know that:
$\gamma\left(k\right) = \mathrm{E}\left[y_t,y_{t-k}\right]$
so I can write:
$\gamma\left(k\right) = \phi_1 \mathrm{E}\left[y_{t-1}y_{t-k}\right]+\phi_2 \mathrm{E}\left[y_{t-2}y_{t-k}\right]+\theta_1 \mathrm{E}\left[\epsilon_{t-1}y_{t-k}\right]+\mathrm{E}\left[\epsilon_{t}y_{t-k}\right]$
then, to derive the analytic version of the autocovariance function, I need to substitute values of $k$ - 0, 1, 2 ... until I get a recursion that is valid for all $k$ greater than some integer.
Therefore, I substitute $k=0$ and work this through to get:
$$
\gamma \left(0\right) = \mathrm{E}\left[y_t,y_t\right]
 =  \phi_1 \mathrm{E}\left[y_{t-1}y_t\right] + \phi_2 \mathrm{E}\left[y_{t-2}y_t\right]+\theta_1 \mathrm{E}\left[\epsilon_{t-1}y_t\right]+\mathrm{E}\left[\epsilon_ty_t\right]\\
$$
now I can simplify the first two of these terms, and then substitute for $y_t$ as before:
$$
\gamma\left(0\right) = \phi_1 \gamma\left(1\right) + \phi_2 \gamma\left(2\right)\\
+ \theta_1 \mathrm{E}\left[\epsilon_{t-1} \left(\phi_1 y_{t-1} +\phi_2 y_{t-2} +\theta_1 \epsilon_{t-1} + \epsilon_t \right)\right]\\
+ \mathrm{E}\left[\epsilon_t \left(\phi_1 y_{t-1} +\phi_2 y_{t-2} +\theta_1 \epsilon_{t-1} + \epsilon_t \right)\right]
$$
then I multiply out the eight terms, which are:
$$
+\theta_1\phi_1\mathrm{E}\left[\epsilon_{t-1}y_{t-1}\right]\\
+\theta_1\phi_2\mathrm{E}\left[\epsilon_{t-1}y_{t-2}\right]\\
+\theta_1^2\mathrm{E}\left[\left(\epsilon_{t-1}\right)^2\right]=\theta_1^2\sigma_{\epsilon}^2\\
+\theta_1\mathrm{E}\left[\epsilon_{t-1}\epsilon_{t}\right]=\theta_1\mathrm{E}\left[\epsilon_{t-1}\right]\mathrm{E}\left[\epsilon_{t}\right]=0\\
+\phi_1\mathrm{E}\left[\epsilon_{t}y_{t-1}\right]\\
+\phi_2\mathrm{E}\left[\epsilon_{t}y_{t-2}\right]\\
+\theta_1\mathrm{E}\left[\epsilon_t\epsilon_{t-1}\right] = \theta_1\mathrm{E}\left[\epsilon_{t}\right]\mathrm{E}\left[\epsilon_{t-1}\right]=0 \\
+\mathrm{E}\left[\left(\epsilon_t\right)^2\right] = \sigma_{\epsilon}^2
$$
So, I am left needing to resolve the four remaining terms.  I want to use the same logic for lines 1, 2, 5 and 6 as I used on lines 4 and 7 - for example for line 1:
$\theta_1\phi_1\mathrm{E}\left[\epsilon_{t-1}y_{t-1}\right] = \theta_1\phi_1\mathrm{E}\left[\epsilon_{t-1}\right]\mathrm{E}\left[y_{t-1}\right] = 0$ because $\mathrm{E}\left[\epsilon_{t-1}\right]=0$.
Similarly for lines 2, 5 and 6.  But I have a model solution that suggests the expression for $\gamma\left(0\right)$ simplifies to:
$\gamma\left(0\right) = \phi_1\gamma\left(1\right)+\phi_2\gamma\left(2\right) +\theta_1\left(\phi_1+\theta_1\right)\sigma_{\epsilon}^2+\sigma_{\epsilon}^2$
This suggests my simplification as described above would miss the term with the coefficient $\phi_1$ - which under my logic should be 0. Is my logic at fault, or is the model solution I found incorrect?
The worked solution also suggest that "analogously" $\gamma\left(1\right)$ can be found as:
$\gamma\left(1\right) = \phi_1\gamma\left(0\right)+\phi_2\gamma\left(1\right) + \theta_1\sigma_{\epsilon}^2$
and for $k>1$:
$\gamma\left(k\right) = \phi_1\gamma\left(k-1\right)+\phi_2\left(k-2\right)$
I hope the question is clear.  Any assistance will be much appreciated. Thank you in advance.
This is a question related to my research, and is not in preparation for any exam or coursework.
 A: If the ARMA process is causal there is a general formula that provides the autocovariance coefficients.
Consider the causal $\text{ARMA}(p,q)$ process
$$
y_t = \sum_{i = 1}^p \phi_i y_{t-1} + \sum_{j = 1}^q \theta_j \epsilon_{t - j} + \epsilon_t,
$$
where $\epsilon_t$ is a white noise with mean zero and variance $\sigma_\epsilon^2$.
By the causality property, the process can be written as
$$
y_t = \sum_{j = 0}^\infty \psi_j \epsilon_{t - j},
$$
where $\psi_j$ denotes the $\psi$-weights.
The general homogeneous equation for the autocovariance coefficients of a causal $\text{ARMA}(p,q)$ process is
$$
\gamma (k) - \phi_1 \gamma (k-1) - \cdots - \phi_p \gamma (k-p) = 0,
\quad k \geq \max (p, q+1),
$$
with initial conditions
$$
\gamma (k) - \sum_{j = 1}^p \phi_j \gamma (k-j)
  = \sigma_\epsilon^2 \sum_{j = k}^q \theta_j \psi_{j - k},
\quad 0 \leq k < \max (p, q+1).
$$
A: Your calculation mistake in your original question lies in 
$$\theta_1\phi_1\mathrm{E}\left[\epsilon_{t-1}y_{t-1}\right] = \theta_1\phi_1\mathrm{E}\left[\epsilon_{t-1}\right]\mathrm{E}\left[y_{t-1}\right] = 0 \qquad \text{(mistaken)}$$
You cannot separate the expectation $\mathrm{E}\left[\epsilon_{t-1}y_{t-1}\right]$ - $\epsilon_{t-1}$ and $y_{t-1}$ are not independent.
A: I would like to share here an approach for calculating the initial conditions mentioned in the answer from @QuantIbex, since the appearance of the $\psi$ parameters from the $MA(\infty)$ representation in that equation does not allow for immediate calculation based on the original ARMA parameters.
Using ARAM(2,2) as an example, we first multiply the ARMA equation
$$\phi_0 y_t + \phi_1 y_{t-1} + \phi_2 y_{t-2} = \theta_0 \epsilon_t + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2}$$
by $\epsilon_t$, $\epsilon_{t-1}$, and $\epsilon_{t-2}$, and take expectations to get three equations which can be written is a matrix form as
$$
\begin{bmatrix}
\delta_0 & 0 & 0\\
 \delta_1 & \delta_0 & 0\\
\delta_2 & \delta_1 & \delta_0
\end{bmatrix}
\begin{bmatrix}
\phi_0 \\ \phi_1 \\ \phi_2
\end{bmatrix}
=
\sigma^2_{\epsilon}
\begin{bmatrix}
\theta_0 \\ \theta_1 \\ \theta_2 
\end{bmatrix}
$$
where $\delta_{j-i} = \mathbb{E}[\epsilon_{t-j}, y_{t-i}]$ and which is equal to 0 for $j > i$. The left-hand side of the above can be re-arranged so that we can write
$$
\mathbf{L} \boldsymbol{\delta} = 
\begin{bmatrix}
\phi_0 & 0 & 0\\
 \phi_1 & \phi_0 & 0\\
\phi_2 & \phi_1 & \phi_0
\end{bmatrix}
\begin{bmatrix}
\delta_0 \\ \delta_1 \\ \delta_2
\end{bmatrix}
=
\sigma^2_{\epsilon}
\begin{bmatrix}
\theta_0 \\ \theta_1 \\ \theta_2 
\end{bmatrix}
=
\sigma^2_{\epsilon}
\boldsymbol{\theta} 
$$
which can be readily solved by forward substitution which we express formally as
$$
\boldsymbol{\delta} = \sigma^2_{\epsilon} \boldsymbol{L}^{-1} \boldsymbol{\theta}.  
$$
Often we would have $\phi_0 = \theta_0 = 1$, which would results in $\delta_0 = \sigma^2_{\epsilon}$.
Then we multiply the ARMA equation by $y_t$, $y_{t-1}$, and $y_{t-2}$, and take expectations to get three equations which can be written is a matrix form as
$$
\mathbf{\Gamma} \boldsymbol{\phi} = 
\begin{bmatrix}
\gamma_0 & \gamma_1 & \gamma_2 \\
\gamma_1 & \gamma_0 & \gamma_1 \\ 
\gamma_2 & \gamma_1 & \gamma_0 
\end{bmatrix}
\begin{bmatrix}
\phi_0 \\ \phi_1 \\ \phi_2
\end{bmatrix}
=
\begin{bmatrix}
\delta_0 & \delta_1 & \delta_2 \\
 0 & \delta_0 & \delta_1 \\
 0 & 0 & \delta_0
\end{bmatrix}
\begin{bmatrix}
\theta_0 \\ \theta_1 \\ \theta_2 
\end{bmatrix}
$$
where
$$\gamma_{|i-j|} = \mathbb{E}[y_{t-i}y_{t-j}] = \mathbb{Cov}[y_{t-i},y_{t-j}]$$
where the last equality is due to us using a "demeaned" process for which $\mathbb{E}[y_{t}] = 0$. The above can also be rearranged as
$$
\mathbf{\Phi} \boldsymbol{\gamma} = 
\begin{bmatrix}
\phi_0 & \phi_1 & \phi_2 \\
\phi_1 & \phi_0 + \phi_2 & 0\\ 
\phi_2 & \phi_1 & \phi_0 
\end{bmatrix}
\begin{bmatrix}
\gamma_0 \\ \gamma_1 \\ \gamma_2
\end{bmatrix}
=
\begin{bmatrix}
\theta_0 & \theta_1 & \theta_2 \\
\theta_1 & \theta_2 & 0 \\
\theta_2 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
\delta_0 \\ \delta_1 \\ \delta_2
\end{bmatrix}
=
\mathbf{U} \boldsymbol{\delta}
$$
which leads to a solution for the initial condition of the auto-covariance function
$$
\boldsymbol{\gamma} = 
\mathbf{\Phi}^{-1}\mathbf{U} \boldsymbol{\delta} =
\sigma^2_{\epsilon} \mathbf{\Phi}^{-1}\mathbf{U} \boldsymbol{L}^{-1} \boldsymbol{\theta}.
$$
The general forms for the triangular matrices $\mathbf{L}$ and $\mathbf{U}$ are apprant from the above example. The elements of $\mathbf{\Phi}$ in row $i$ and column $j$ can be inferred by looking at the row numbers of $\gamma_i$ in column $j$ of the $\boldsymbol{\Gamma}$ matrix and are in general given by (I believe)
$$
(\mathbf{\Phi})_{ij} = \sum_{k=0}^{p} \phi_k \mathbb{I}(|k - i| = j)
$$
where we assume that the row and column index begins with $0$ (not $1$) and where $\mathbb{I}(A)$ is the indicator function which is equal to 1 when $A$ is TRUE and $0$ otherwise.
A: OK. So the process of writing the post actually pointed me to the solution.
Consider the Expectation terms 1, 2, 5 and 6 from above that I thought should be 0.  
Immediately for terms 5 - $\mathrm{E}\left[\epsilon_ty_{t-1}\right]$ - and 6 - $\mathrm{E}\left[\epsilon_ty_{t-2}\right]$:  these terms are definitely zero, because $y_{t-1}$ and $y_{t-2}$ are independent of $\epsilon_t$ and $\mathrm{E}\left[\epsilon_t\right] = 0$.
However, terms 1 and 2 look as though the Expectation is of two correlated variables.  So, consider the expressions for $y_{t-1}$ and $y_{t-2}$ thus:
$$
y_{t-1} = \phi_1y_{t-2}+\phi_2y_{t-3}+\theta_1\epsilon_{t-2}+\epsilon_{t-1}\\
y_{t-2} = \phi_1y_{t-3}+\phi_2y_{t-4}+\theta_1\epsilon_{t-3}+\epsilon_{t-2}
$$
And recall term 1 - $\phi_1\theta_1\mathrm{E}\left[\epsilon_{t-1}y_{t-1}\right]$. If we multiply both sides of the expression for $y_{t-1}$ by $\epsilon_{t-1}$ and then take Expectations, it is clear that all terms on the right hand side except the last become zero (because the values of $y_{t-2}$, $y_{t-3}$, and $\epsilon_{t-2}$ are independent of $\epsilon_{t-1}$ and $\mathrm{E}\left[\epsilon_{t-1}\right]=0$) to give:
$$\mathrm{E}\left[\epsilon_{t-1}y_{t-1}\right] = \mathrm{E}\left[\left(\epsilon_{t-1}\right)^2\right] = \sigma_{\epsilon}^2$$    
So term 1 becomes $+\phi_1\theta_1\sigma_{\epsilon}^2$. For term 2, it should be clear that, by the same logic, all terms are zero.
Hence the original model answer was correct.
However, if anyone can suggest an alternative way to obtain a general (even if messy) solution, I would be very pleased to hear it!
