convergence in probability of mean of sequence $X_{k+1}=\beta X_k+\epsilon_k$ Let $X_1=0$ and $X_{k+1}=\beta X_k+\epsilon_k$ with $\epsilon_k$ iid normally distributed, $\epsilon_k\sim N(\mu,\sigma^2)$. Let $\mu$ and $\sigma^2$ be fixed.
How can you show the mean of $X_k$ for $k=1,...,n$ converges, e.g. $$\frac1n\sum_{k=1}^nX_k\xrightarrow{P}{} \frac{\mu}{1-\beta} \;?$$
 A: The final version is based on the helpful (and constructive) comments (see below) that pointed out mistakes and flaws of the previous efforts to answer the question.

By substituting iteratively you find that $$\begin{align*}X_{k+1}&=βX_{k}+\epsilon_k=β(βX_{k-1}+\epsilon_{k-1})+\epsilon_k=β^2X_{k-1}+β\epsilon_{k-1}+\epsilon_k=\ldots\\\\&=β^kX_1+\left(\sum_{l=0}^{k-1}β^l\epsilon_{k-l}\right)\\&=\sum_{l=0}^{k-1}β^l\epsilon_{k-l}\end{align*}$$ since $X_1=0$. Thus \begin{array}{rcrcrr}X_2&=&\epsilon_1\\X_3&=&β\epsilon_1&+&\epsilon_2\\X_4&=&β^2\epsilon_1&+&β\epsilon_2&+&\epsilon_3&\\\ldots&&\ldots&&\ldots&&\ldots\\X_n+1&=&β^{n-1}\epsilon_1&+&β^{n-2}\epsilon_2&+&β^{n-3}\epsilon_3&+\ldots+β\epsilon_{n-1}+\epsilon_n\end{array}
Summing up both sides, we obtain $$\begin{align*}\sum_{k=1}^{n}X_{k+1}&=\frac{1-β^n}{1-β}\epsilon_1+\frac{1-β^{n-1}}{1-β}\epsilon_2+\ldots+\frac{1-β^2}{1-β}\epsilon_{n-1}+\frac{1-β}{1-β}\epsilon_n\\&=\sum_{k=1}^{n}\frac{1-β^{n+1-k}}{1-β}\epsilon_k=\ldots\\&=\frac{1}{1-β}\left(\sum_{k=1}^{n}\epsilon_k-β^{n+1}\sum_{k=1}^{n}\frac{1}{β^k}\epsilon_k\right)\end{align*}$$ and finally $$\bar{X}_n=\frac{1}{n(1-β)}\left(\sum_{k=1}^{n}\epsilon_k-β^{n+1}\sum_{k=1}^{n}\frac{1}{β^k}\epsilon_k\right)$$ From this form we can calculate $E[\bar{X}_n]$ and (especially) $Var(\bar{X}_n)$ since the $\epsilon_k$ are independent. Indeed $$E[\bar{X}_n]=\frac{1}{n(1-β)}\left(\sum_{k=1}^{n}μ-β^{n+1}\sum_{k=1}^{n}\frac{1}{β^k}μ\right)=\frac{μ}{1-β}\left(1-\frac{β^{n+1}-1}{n(1-β)}\right) \longrightarrow \frac{μ}{1-β}$$ as $n \to \infty$ and similarly $$\begin{align*}Var[\bar{X}_n]&=\frac{1}{n^2(1-β)^2}\left(\sum_{k=1}^{n}σ^2-(β^2)^{n+1}\sum_{k=1}^{n}\frac{1}{(β^2)k}σ^2\right)\\\\&=\frac{σ^2}{(1-β)^2}\left(\frac{1}{n}-\frac{(β^2)^{n+1}-1}{n(1-β^2)}\right) \, \longrightarrow \, 0\end{align*}$$ as $n \to \infty$. Hence $$E[\bar{X}_n^2]=Var(\bar{X}_n)+E[\bar{X}_n]^2 \longrightarrow 0+\frac{μ^2}{(1-β)^2}$$ which is enough to obtain $\mathcal L^2$ convergence $$E\left[\left(X_n-\frac{μ}{1-β}\right)^2\right]=E\left[\bar{X}_n^2-2\bar{X}_n\frac{μ}{1-β}+\frac{μ^2}{(1-β)^2}\right]\to \frac{2μ^2}{(1-β)^2}-\frac{2μ^2}{(1-β)^2}=0$$ from which you can conclude that $$\overline{X}_n=\frac{1}{n}\sum_{k=1}^{n}X_k\overset{p}\rightarrow \frac{μ}{1-β}$$ since convergence in $\mathcal L^r$ for $1\le r$ implies convergence in probability (see here).

Summing up both sides of the given relation $$\sum_{k=1}^{n}X_{k+1}=β\sum_{k=1}^{n}X_k+\sum_{k=1}^{n}\epsilon_k$$ gives you after simple calculations that $$\bar{X}_n=\frac{1}{n(1-β)}\left(\sum_{k=1}^{n}\epsilon_k-X_{n+1}\right)$$ and thus obtaining a closed form for $X_{n+1}$ (as above) is enough to obtain the result with less calculations than above.
