Definitions
The ARMA model
$$x_n=-\sum_{p=1}^P a_px_{n-p}+\sum_{q=0}^Qb_qw_{n-q} \tag{1}$$
where $w_n$ is zero mean stationary white noise with unit variance.
Question
To derive the autocorrelation function of the ARMA process, $R_{XX}[r]=E[x_nx_{n+r}]$ I think we can do this two ways:
- Substitute eqn 1 into $x_{n+r}$
$$R_{XX}[r]=E\left[x_n\left\{-\sum_{p=1}^Pa_px_{n+r-p}+\sum_{q=0}^Qb_qw_{n+r-q}\right\}\right]$$
$$=-\sum_{p=1}^Pa_pE[x_nx_{n+r-p}]+\sum_{q=0}^Qb_qE[x_nw_{n+r-q}]$$
$$=-\sum_{p=1}^Pa_pR_{XX}[r-p]+\sum_{q=0}^Qb_qR_{XW}[r-q]$$
- Substitute eqn 1 into $x_{n}$
$$R_{XX}[r]=E\left[\left\{-\sum_{p=1}^Pa_px_{n-p}+\sum_{q=0}^Qb_qw_{n-q}\right\}x_{n+r}\right]$$
$$=-\sum_{p=1}^Pa_pE[x_{n-p}x_{n+r}]+\sum_{q=0}^Qb_qE[w_{n-q}x_{n+r}]$$
$$=-\sum_{p=1}^Pa_pR_{XX}[r+p]+\sum_{q=0}^Qb_qR_{XW}[-r-q]$$
But $R_{XX}[r-p]\ne R_{XX}[r+p]$ and $R_{XW}[r-q]\ne R_{XX}[r+q]$, so what mistake have I made? Is only one method correct, and why?
Both ways are correct, but only one way is useful.
Before, be careful with your definition of cross-correlation ($R_{XW}[t]\ne R_{WX}[t]$), it's a matter of convention, but you must be consistent. If we define $R_{XW}[t]=E[x_n w_{n-t}]$, then the first derivation has an error, it should be: $$R_{XX}[r]=-\sum_{p=1}^Pa_pR_{XX}[r-p]+\sum_{q=0}^Qb_qR_{XW}[q-r]\tag{1}$$ The second is ok: $$R_{XX}[r]=-\sum_{p=1}^Pa_pR_{XX}[r+p]+\sum_{q=0}^Qb_qR_{XW}[r+q]\tag{2}$$
Both are correct (and the equivalence is non-trivial) but only $(1)$ is useful, in the Yule-Walker context (we want to write $R_{XX}[0]$, $R_{XX}[1]$ $\cdots$ and "solve" the system), because only the first one allows us to exploit causality, the fact that $x_n$ does not depend (is independent of) future values of $w$, that is, $E[x_n w_{n-t}]=R_{XW}[t]=0 $ for $t<0$.
Once more, I encourage you to get familiar with the AR case, which is conceptually similar, but simpler.