Let $Z_n \sim \mathcal{N}(0,1)$ and $$ X_{n+1} = \frac{1}{2} X_n - \frac{5}{16} X_{n-1} + Z_n. $$
Find the limiting joint distribution of $X_n$ and $X_{n-1}$ as $n \to \infty$.
Attempt
My approach was to rewrite the system as the one-lag system \begin{align*} \begin{pmatrix} X_{n+1} \\ X_n \end{pmatrix} &= \begin{pmatrix} \frac{1}{2} & -\frac{5}{16} \\ 1 & 0 \end{pmatrix} \begin{bmatrix} X_n \\ X_{n-1} \end{bmatrix} + \begin{pmatrix} 1 \\ 0 \end{pmatrix} Z_n \\ &=: A \begin{bmatrix} X_n \\ X_{n-1} \end{bmatrix} + B Z_n \end{align*} and solve for the covariance using the formula $$ C = \sum_{k=0}^\infty A^k BB^t (A^k)^t. $$ However, I don't know if this is the most reasonable way to solve for the covariance. It's quite easy to plug into MATLAB and evaluate the series but I wonder how to solve for $C$ without resorting to numerical methods.