[proofplan]
We prove the recursion by writing the $k$-variable prediction error as a linear combination of two $(k-1)$-variable errors: the forward error obtained by predicting $X_{k+1}$ from $X_k,\dots,X_2$, and the backward error obtained by predicting $X_1$ from $X_2,\dots,X_k$. Stationarity implies that the reversed Durbin-Levinson coefficients give the backward predictor and that both forward and backward error variances equal $v_{k-1}$. We choose the scalar multiplying the backward error so that the resulting error is orthogonal to $X_1$, and this gives the formula for $\phi_{k,k}$. Expanding the resulting error gives the coefficient recursion, and computing its variance gives $v_k = v_{k-1}(1-\phi_{k,k}^2)$.
[/proofplan]
[step:Write the normal equations for the forward prediction problem]
For $k \in \mathbb{N}$, define the forward prediction error
\begin{align*}
E_k
=
X_{k+1}-\sum_{j=1}^k \phi_{k,j}X_{k+1-j}.
\end{align*}
Since $\Gamma_k$ is nonsingular, the least-squares coefficients are unique. The orthogonality condition for best linear prediction is
\begin{align*}
\mathbb{E}[E_k X_{k+1-i}] = 0,
\qquad 1 \le i \le k.
\end{align*}
Using stationarity, this is equivalent to
\begin{align*}
\gamma(i)-\sum_{j=1}^k \phi_{k,j}\gamma(i-j)=0,
\qquad 1 \le i \le k.
\end{align*}
For $k=1$, the equation gives
\begin{align*}
\phi_{1,1}=\frac{\gamma(1)}{\gamma(0)}=\frac{\gamma(1)}{v_0},
\end{align*}
and the variance identity follows from a direct calculation:
\begin{align*}
v_1
&=
\mathbb{E}\left[\left(X_2-\phi_{1,1}X_1\right)^2\right] \\
&=
\gamma(0)-2\phi_{1,1}\gamma(1)+\phi_{1,1}^2\gamma(0) \\
&=
v_0(1-\phi_{1,1}^2).
\end{align*}
Hence assume $k \ge 2$ for the remaining argument.
[/step]
[step:Construct the forward and backward errors of order $k-1$]
Define the shifted forward error
\begin{align*}
A_k
=
X_{k+1}-\sum_{j=1}^{k-1}\phi_{k-1,j}X_{k+1-j}.
\end{align*}
This is the prediction error for $X_{k+1}$ from $\operatorname{span}\{X_k,\dots,X_2\}$, obtained from the order $k-1$ predictor by shifting time forward by one. Therefore
\begin{align*}
\mathbb{E}[A_k X_{k+1-i}]=0,
\qquad 1 \le i \le k-1,
\end{align*}
and
\begin{align*}
\mathbb{E}[A_k^2]=v_{k-1}.
\end{align*}
Define the backward error
\begin{align*}
B_k
=
X_1-\sum_{j=1}^{k-1}\phi_{k-1,k-j}X_{k+1-j}.
\end{align*}
We verify that $B_k$ is orthogonal to $\operatorname{span}\{X_k,\dots,X_2\}$. For $1 \le i \le k-1$,
\begin{align*}
\mathbb{E}[B_k X_{k+1-i}]
&=
\gamma(k-i)-\sum_{j=1}^{k-1}\phi_{k-1,k-j}\gamma(i-j).
\end{align*}
Set $r=k-i$. Then $1 \le r \le k-1$, and after replacing $j$ by $k-\ell$,
\begin{align*}
\mathbb{E}[B_k X_{k+1-i}]
&=
\gamma(r)-\sum_{\ell=1}^{k-1}\phi_{k-1,\ell}\gamma(r-\ell).
\end{align*}
This is the $r$-th normal equation for the order $k-1$ predictor, hence it is zero. Thus $B_k$ is the backward prediction error for $X_1$ from $\operatorname{span}\{X_k,\dots,X_2\}$, and by stationarity its variance is
\begin{align*}
\mathbb{E}[B_k^2]=v_{k-1}.
\end{align*}
[guided]
The predictor of order $k-1$ can be shifted because the process is stationary. The usual order $k-1$ error predicts $X_k$ from $X_{k-1},\dots,X_1$. Shifting every index up by one gives
\begin{align*}
A_k
=
X_{k+1}-\sum_{j=1}^{k-1}\phi_{k-1,j}X_{k+1-j},
\end{align*}
which predicts $X_{k+1}$ from $X_k,\dots,X_2$. Since covariances depend only on time differences, the same normal equations hold, so $A_k$ is orthogonal to each of $X_k,\dots,X_2$, and its variance is $v_{k-1}$.
The backward error uses the same coefficients in reverse order:
\begin{align*}
B_k
=
X_1-\sum_{j=1}^{k-1}\phi_{k-1,k-j}X_{k+1-j}.
\end{align*}
We check the orthogonality directly. For $1 \le i \le k-1$,
\begin{align*}
\mathbb{E}[B_k X_{k+1-i}]
&=
\mathbb{E}[X_1X_{k+1-i}]
-
\sum_{j=1}^{k-1}\phi_{k-1,k-j}\mathbb{E}[X_{k+1-j}X_{k+1-i}] \\
&=
\gamma(k-i)-\sum_{j=1}^{k-1}\phi_{k-1,k-j}\gamma(i-j).
\end{align*}
Put $r=k-i$. Reindex the sum by $\ell=k-j$. Then
\begin{align*}
\mathbb{E}[B_k X_{k+1-i}]
&=
\gamma(r)-\sum_{\ell=1}^{k-1}\phi_{k-1,\ell}\gamma(r-\ell).
\end{align*}
This is exactly the $r$-th normal equation for the order $k-1$ predictor, so it equals $0$. Thus $B_k$ is orthogonal to $X_k,\dots,X_2$. Stationarity also gives $\mathbb{E}[B_k^2]=v_{k-1}$, because the joint covariance matrix of $(X_1,X_2,\dots,X_k)$ is the same Toeplitz matrix as the reversed prediction problem of order $k-1$.
[/guided]
[/step]
[step:Choose the reflection coefficient by enforcing orthogonality to $X_1$]
Since $A_k$ is orthogonal to $X_k,\dots,X_2$ and $B_k$ is orthogonal to $X_k,\dots,X_2$, every random variable of the form $A_k-\alpha B_k$ with $\alpha \in \mathbb{R}$ is orthogonal to $X_k,\dots,X_2$. We choose $\alpha$ so that it is also orthogonal to $X_1$.
First compute
\begin{align*}
\mathbb{E}[A_k X_1]
&=
\gamma(k)-\sum_{j=1}^{k-1}\phi_{k-1,j}\gamma(k-j).
\end{align*}
Since $B_k$ is orthogonal to $\operatorname{span}\{X_k,\dots,X_2\}$ and
\begin{align*}
X_1
=
B_k+\sum_{j=1}^{k-1}\phi_{k-1,k-j}X_{k+1-j},
\end{align*}
we have
\begin{align*}
\mathbb{E}[B_kX_1]=\mathbb{E}[B_k^2]=v_{k-1}.
\end{align*}
Thus the condition
\begin{align*}
\mathbb{E}[(A_k-\alpha B_k)X_1]=0
\end{align*}
is equivalent to
\begin{align*}
\alpha
=
\frac{\mathbb{E}[A_kX_1]}{v_{k-1}}
=
\frac{\gamma(k)-\sum_{j=1}^{k-1}\phi_{k-1,j}\gamma(k-j)}{v_{k-1}}.
\end{align*}
Define
\begin{align*}
\phi_{k,k}
=
\frac{\gamma(k)-\sum_{j=1}^{k-1}\phi_{k-1,j}\gamma(k-j)}{v_{k-1}}.
\end{align*}
Then $A_k-\phi_{k,k}B_k$ is orthogonal to all of $X_k,\dots,X_1$.
[/step]
[step:Expand the new error to identify the coefficient recursion]
Expanding $A_k-\phi_{k,k}B_k$ gives
\begin{align*}
A_k-\phi_{k,k}B_k
&=
X_{k+1}
-\sum_{j=1}^{k-1}\phi_{k-1,j}X_{k+1-j}
-\phi_{k,k}X_1
+\phi_{k,k}\sum_{j=1}^{k-1}\phi_{k-1,k-j}X_{k+1-j} \\
&=
X_{k+1}
-\sum_{j=1}^{k-1}\left(\phi_{k-1,j}-\phi_{k,k}\phi_{k-1,k-j}\right)X_{k+1-j}
-\phi_{k,k}X_1.
\end{align*}
The random variable on the left is orthogonal to $\operatorname{span}\{X_k,\dots,X_1\}$. By nonsingularity of $\Gamma_k$, the best linear predictor from this span is unique, so
\begin{align*}
E_k=A_k-\phi_{k,k}B_k.
\end{align*}
Comparing coefficients of $X_{k+1-j}$ for $1 \le j \le k-1$ and of $X_1=X_{k+1-k}$ gives
\begin{align*}
\phi_{k,j}
=
\phi_{k-1,j}-\phi_{k,k}\phi_{k-1,k-j},
\qquad 1 \le j \le k-1.
\end{align*}
Together with the definition of $\phi_{k,k}$, this proves the coefficient recursion.
[/step]
[step:Compute the variance of the updated prediction error]
Since $E_k=A_k-\phi_{k,k}B_k$, we compute
\begin{align*}
v_k
&=
\mathbb{E}[E_k^2] \\
&=
\mathbb{E}[A_k^2]
-2\phi_{k,k}\mathbb{E}[A_kB_k]
+\phi_{k,k}^2\mathbb{E}[B_k^2].
\end{align*}
Because $A_k$ is orthogonal to $X_k,\dots,X_2$ and
\begin{align*}
B_k
=
X_1-\sum_{j=1}^{k-1}\phi_{k-1,k-j}X_{k+1-j},
\end{align*}
we have
\begin{align*}
\mathbb{E}[A_kB_k]=\mathbb{E}[A_kX_1].
\end{align*}
By the definition of $\phi_{k,k}$,
\begin{align*}
\mathbb{E}[A_kX_1]=\phi_{k,k}v_{k-1}.
\end{align*}
Using also $\mathbb{E}[A_k^2]=\mathbb{E}[B_k^2]=v_{k-1}$, we obtain
\begin{align*}
v_k
&=
v_{k-1}
-2\phi_{k,k}^2v_{k-1}
+\phi_{k,k}^2v_{k-1} \\
&=
v_{k-1}(1-\phi_{k,k}^2).
\end{align*}
This is the stated variance recursion, and the proof is complete.
[/step]
