Ok, so the plan of the proof is as follows:
First, I am going to fix the Soong and Juang's proof that if all roots of $A$ lie inside the unit circle, then all roots of $P$ and $Q$ lie on the unit circle. (I am going to follow their notation, too.)
Then, I am going to prove that all roots of the LPC polynomial are indeed inside the unit circle.
Finally, I'll prove that the roots of $P$ and $Q$ are intertwined.
The first part is easy. We define
$$
H(z) \triangleq z^{-(M+1)}\frac{A(z^{-1})}{A(z)}
$$
Since A is a polynomial with real coefficients, complex conjugation is a permutation on the set of its roots, therefore H can be factored as
$$
H(z) = z^{-1} \prod_{i=1}^M {\frac{z_iz-1}{z-\overline{z_i}}},
$$
where $z_i$ are the roots of $A_M(z)$.
Now, since
$ |z_iz - 1|^2 - |z-\overline{z_i}|^2 = (z_iz - 1)(\overline{z_iz} - 1) - (z - \overline{z_i})(\overline{z} - z_i) = |z_i|^2|z|^2 + 1 - zz_i -\overline{zz_i} - |z|^2 - |z_i|^2 + zz_i + \overline{z}\overline{z_i} = (1-|z|^2)(1-|z_i|^2)$,
Therefore,
$$
|H(z)| \text{ is } \begin{cases}
>1 &\quad\text{if } |z| < 1\\
<1 &\quad\text{if } |z| > 1\\
=1 &\quad\text{if } |z| = 1\\
\end{cases}
$$
Because the solution to P(z) = 0 or Q(z) = 0 requires that H(z) = ± 1, we conclude that P(z) and Q(z) can only have zeros on the unit circle.
OK, now how do we prove that all zeros of $A(z)$ lie inside the unit circle? We are going to prove that the polynomials $A_M(z)$ satisfy a recursive relation:
$$
A_{M+1}(z) = A_M(z) - \kappa_M z^{-(M+1)}A_M(z^{-1}),
$$
where $|\kappa_M| < 1$. ($\kappa_M$ are called reflection coefficients.) Then we'll prove recursively that all roots of $A_M(z)$ are in the unit circle.
$A_M(z) = 1 + a_{M,1} z^{-1} + \ldots + a_{M,M} z^{-M}$ is a polynomial whose coefficients minimize the forward prediction error:
$$
E_{f,M} = \sum_k{\left(x(k) + \sum_{m=1}^M {x(k-m) a_{M,m}}\right)^2}
$$
If we write down the optimality condition for $a_{M,m}$, we get the Wiener-Hopf equations:
$$
R_M \vec{a_M} = - \vec{r_{f,M}},
$$
where
$$
\vec{a_M} \triangleq [a_1, a_2, \ldots, a_M]^T,\\
\vec{r_{f,M}} \triangleq [r(1), r(2), \ldots, r(M)]^T,\\
r(l) \triangleq \sum_k{x(k) x(k-l)}, \\
R_M \triangleq \begin{pmatrix}
r(0) & r(1) & \ldots & r(M-1) \\
r(1) & r(0) & \ldots & r(M-2) \\
\vdots & \vdots & \ddots & \vdots \\
r(M-1) & r(M-2) & \ldots & r(0) \\
\end{pmatrix}
$$
Similarly, the solution to the problem of backwards linear prediction
$$
\text{minimize }
E_{b,M} = \sum_k{\left(x(k-M) + \sum_{m=1}^M {x(k-m+1) b_{M,m}}\right)^2}
$$
must satisfy the equation
$$
R_M \vec{b_M} = - \vec{r_{b,M}},
$$
where
$$
\vec{r_{b,M}} \triangleq [r(M), r(M-1), \ldots, r(1)]^T
$$
Define the coidentity matrix $J_M$ as
$$
J_M = \begin{pmatrix}
0 & 0 & \ldots & 0 & 1 \\
0 & 0 & \ldots & 1 & 0 \\
\vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 1 & \ldots & 0 & 0 \\
1 & 0 & \ldots & 0 & 0 \\
\end{pmatrix}
$$
It is easy to see that
$$
J_MR_M = R_MJ_M,
$$
therefore
$$
R_M \vec{a_M} = -\vec{r_{f,M}} = -J_M\vec{r_{b,M}} = J_M R_M \vec{b_M} =
R_M J_M \vec{b_M}.
$$
Since, as I have shown, the matrix $R_M$ is nonsingular, it follows that $\vec{a_M} = J_M \vec{b_M}$.
Furthermore,
$$
E_{f,M} = r(0) + \vec{r_{f,M}^T} \vec{a_M} = r(0) + \vec{r_{f,M}^T}J_M J_M \vec{a_M} = r(0) + \vec{r_{b,M}^T} \vec{b_M} = E_{b,M} = E_M
$$
Consider the following expression:
$$
\begin{pmatrix}
R_M & \vec{r_{b,M}} \\
\vec{r_{b,M}}^T & r(0) \\
\end{pmatrix}
\begin{pmatrix}
1 \\
\vec{a_{M-1}} \\
0\\
\end{pmatrix} =
\begin{pmatrix}
E_{M-1}\\
0_{(M-1)\times 1} \\
K_M\\
\end{pmatrix}
$$
where
$$
K_M = r(M) + \vec{r_{b,M-1}^T}\vec{a_{M-1}} = r(M) + \vec{r_{b,M-1}^T}J_M\vec{b_{M-1}} = r(M) + \vec{r_{f,M-1}^T}\vec{b_{M-1}}
$$
We define the reflection coefficient as
$$
\kappa_M = \frac{K_M}{E_{M-1}}
$$
From backward linear prediction, we have:
$$
\begin{pmatrix}
r(0) & \vec{r_{f,M}^T} \\
\vec{r_{f,M}} & R_M \\
\end{pmatrix}
\begin{pmatrix}
0 \\
\vec{b_{M-1}} \\
1\\
\end{pmatrix} =
\begin{pmatrix}
K_M\\
0_{(M-1)\times 1} \\
E_{M-1}\\
\end{pmatrix}
$$
Multiplying both sides of this equation by $\kappa_M$ and subtracting it from the previous matrix equation, we get:
$$
R_{M+1}\begin{pmatrix}
1\\
\vec{a_{M-1}} - \kappa_M \vec{b_{M-1}} \\
-\kappa_M\\
\end{pmatrix} =
\begin{pmatrix}
E_{M-1} (1 - \kappa_M^2) \\
0_{M \times 1} \\
\end{pmatrix}
$$
From which it follows that
$$
\vec{a_M} = \begin{pmatrix}
\vec{a_{M-1}} - \kappa_M \vec{b_{M-1}} \\
-\kappa_M\\
\end{pmatrix}
$$
and $E_M = E_{M-1}(1-\kappa_M^2)$. Since $E_M > 0$ and $E_{M-1} > 0$, we get $|\kappa_M| < 1$.
Now let us prove by induction that all roots of $A_M$ lie inside the unit circle.
Base: For $A_1$, $\kappa_1$ is the only root.
Induction step: Suppose we proved that all roots of $A_M$ lie in the unit circle.
$$
A_{M+1}(z) = A_M(z) (1 - \kappa_M H(z)),
$$
So, if $A_{M+1}(z) = 0$, then either $A_M(z) = 0$ or $|H(z)| = 1/\kappa_M > 1$, therefore $|z| < 1$.
The last thing to prove: zeros of P and Q are intertwined. Zeros of P are solutions to $H(z) = -1$ and zeros of Q are solutions to $H(z) = 1$, so it is sufficient to prove that the phase of H monotonously decreases as $z = e^{jw}$ travels the unit circle.
H(z) can be factored as
$$
H(z) = z^{-(M+1)}\prod_{i=1}^M {\frac{1 - \overline{z_i^{-1}}z^{-1}}{1 - z_iz^{-1}}}
$$
We are going to prove that the phase of each multiplier is a monotonously decreasing function of w. An easy way to do this is by using the law of tangents.
$$
\frac{\tan{(\beta - \alpha)/2}}{\tan{(\beta + \alpha)/2}} = \frac{1/r - 1}{1/r + 1}
$$
$$
\angle \frac{1 - \overline{z_i^{-1}}z^{-1}}{1 - z_iz^{-1}} = - 2 \arctan{(\beta - \alpha)/2} = - 2 \arctan{\left(\frac{1/r - 1}{1/r + 1} \tan (90° - (w_i - w)/2)\right)}
$$
is a monotonously decreasing function of w, q.e.d.
