Yule Walker equations of an ARMA(1,1)-process I want to find Yule-Walker equations for a causal ARMA(1,1)-process 
$Y_t=\alpha Y_{t-1}+\epsilon_t+\theta\epsilon_{t-1}$. 
My idea initially was to first multiply both sides by $Y_{t-1}$, so I would get 
$\gamma_0=\alpha\gamma_1+\epsilon_t Y_t+\theta\epsilon_{t-1}Y_t$, where $\gamma_0$ is the variance and $\gamma_1$ is the first autocovariance.
But, I don't know how to continue from there. 
Could someone help me out? Thanks in advance!
 A: Yule Walker (for parameter estimation) is usually only used for AR models, but this method you're using is still a valid technique for finding the autocovariance function. I'm assuming that's what you're after.
Multiply both sides of an ARMA(p,q) model by $Y_{t-k} = \sum_{j=0}^{\infty}\psi_j \epsilon_{t-j-k}$ (this is where causality assumption comes in)
$$
Y_t Y_{t-k} - \phi_1 Y_{t-1}Y_{t-k} - \cdots - \phi_p Y_{t-p}Y_{t-k} = \left[\sum_{l=0}^q \theta_l \epsilon_{t-l}\right]\left[\sum_{j=0}^{\infty}\psi_j \epsilon_{t-j-k} \right].  \tag{1}
$$
Then take expectations. Note that $\psi_j = 0$ for $j < 0$ and $\theta_j = 0$ for $j \not \in \{1,\ldots,q\}$. If $0 \le k < \text{max}(p,q+1)$ you get these $p$ (1 in your case) equations here
$$
\gamma(k) - \phi_1 \gamma(k-1) - \cdots - \phi_p\gamma(k-p) = \sigma^2 \sum_{j=0}^{\infty}\theta_{k+j}\psi_j. \tag{2}
$$
And if $k \ge \text{max}(p,q+1)$ you get
$$
\gamma(k) - \phi_1 \gamma(k-1) - \cdots - \phi_p\gamma(k-p) = 0 \tag{3}
$$
($0$ on RHS because $l > q$).
You can solve your two equations by using results of homogeneous difference equations, or by just solving the system directly.
For your specific ARMA(1,1) model (I'll turn the $\phi$s back into $\alpha$s) you have


*

*$\gamma(0) - \alpha \gamma(1) = \sigma^2[1 + \theta_1 (\theta + \alpha)]$

*$\gamma(1) - \alpha \gamma(0) = \sigma^2[\theta]$

*$\gamma(k) - \alpha \gamma(k-1) = 0$ for $k \ge 2$


Plugging (2) into (1) gives us
$$
\gamma(0) = \sigma^2\left[ 1 + \frac{(\theta+\alpha)^2 }{1 - \alpha^2 } \right]
$$
(which is the same as before) and plugging this into (2) gives us
$$
\gamma(1) = \sigma^2\left[ \theta + \alpha + \frac{(\theta+\alpha)^2\alpha }{1 - \alpha^2 } \right].
$$
If you want to get it at higher lags, just use (3) repeatedly.
When the model is a pure AR model, the RHS of (1) is simplified tremendously, and the resulting set of equations to solve is linear in the parameters $\phi_1,\ldots,\phi_p$.
