Definitions
For an 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.
It is straightforward to show that the autocorrelation of $x_n$, $R_{XX}[r]$ is given by:
$$R_{XX}[r]=-\sum_{p=1}^P a_pR_{XX}[r-p]+c_r \tag{2}$$
where
$$c_r = \begin{cases}\sum_{q=r}^Qb_qh_{q-r} & r\leq Q \\ 0 & r>Q \end{cases} \tag{3}$$
where $h_n$ is the impulse response of the system, i.e. $x_n=\sum_{m=-\infty}^{\infty}h_mw_{n-m}=h_n*w_n$
Question
I have seen in textbooks that the equation for $R_{XX}[r]$ expands to
$$R_{XX}[0]+\sum_{p=1}^Pa_pR_{XX}[-p]=c_0 \tag{4}$$ $$R_{XX}[1]+\sum_{p=1}^Pa_pR_{XX}[1-p]=c_1$$ $$\vdots$$ $$R_{XX}[Q]+\sum_{p=1}^Pa_pR_{XX}[Q-p]=c_Q$$ $$R_{XX}[Q+1]+\sum_{p=1}^Pa_pR_{XX}[Q+1-p]=c_{Q+1}=0$$ $$\vdots$$ $$R_{XX}[Q+P]+\sum_{p=1}^Pa_pR_{XX}[Q+P-p]=0$$
These can then easily be put into the classic matrix form $\mathbf{R}_X'\mathbf{a}'=\mathbf{c}'$.
This makes perfect sense except for one detail - why do the textbooks stop the equations at $R_{XX}[Q+P]$? Why do they not continue to $R_{XX}[Q+P+1]$ etc.?
I have tried fruitlessly to show that $\sum_{p=1}^Pa_pR_{XX}[Q+P+1-p]=0$ meaning $R_{XX}[Q+P+1]=0$ so that there isn't any point in considering these equations.
I also tried to show from $x_n=-\sum_{p=1}^P a_px_{n-p}+\sum_{q=0}^Qb_qw_{n-q}$ that $R_{XX}[Q+P+1]=0$. I would really appreciate it if someone helped me out here, I'm sure the answer is simple and obvious but I just can't see it.
Actually, they don't. Equations (2) and (3) are valid for arbitrarily big lags.
What is true is that the Yule Walker set of equations, when writen in matricial form, "stops" at lag $Q+P$ so that we have $Q+P+1$ rows or equations, which equals the "degrees of freedom" of our ARMA model.
The YW equations, then, should be regarded as a transformation between two sets of parameters that can be used to characterize the model: the filter coefficientes and the autocorrelation values: $$\{a_1,a_2 \cdots a_P, b_0, b_1 \cdots b_Q\} \longleftrightarrow \{ r_0, r_1 \cdots r_{Q+P}\}$$
They provide an answer to two different questions: "Given this ARMA model, which is the autocorrelation?" and "Given this autocorrelation, which is the ARMA model?" - more precisely: which equations should the $(P+Q+1)$ ARMA coefficients verify so that the autocorrelation fits my desired values from lag $0$ to $P+Q$? The point, then, is that you can expect to adjust a model with $(P+Q+1)$ degrees of freedom so that to fit $(P+Q+1)$ values of the autocorrelation.
Let's suppose you have some autocorrelation function $R_x[n]$ in full (or, say, for $n\le 100$), and you wish to find a $(2,4)$ ARMA model that fits it. Then, you have $7$ parameters to find, and you'd write the corresponding YW equations, which only involve $R_x[0]$ to $R_x[6]$. You have then $7$ equations with $7$ unkowns. Assumming you have found one solution (thing are a bit complicated for ARMA, the equations are really non linear), then your ARMA model will generate a process that fits your desired correlation, but only up to $R_x[6]$; for bigger lags, you'd get an autocorrelation given by equations $2,3$ - and this, of course, will not coincide with your own original correlation, in general (unless, perhaps, you knew in advance that your $R_x[n]$ came from an $(2,4)$ ARMA process).
If you are insterested only in the first question (compute autocorrelation given the ARMA model), then indeed you don't need to write the YW equations, you'd just use the eqs. $(2)$, $(3)$ in full.
This is more elegant/easy to grasp in the AR model, in which the relation is one-to-one.