[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*}[/step]

[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: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]