I'm trying to find a reference (including the full formula) for the following. If $X_n = a_1 X_{n-1} + \cdots a_p X_{n-p} + e(n)$ where $\{e(n)\}$ is a white noise, then
$$ X_n=g(e_0,e_1,\ldots,e_n)+\sum_{k=1}^{p}r_k^n\cdot Q_k,\ \mathrm{with}\ \begin{pmatrix} Q_1 \\ Q_2 \\ \vdots \\ Q_p \end{pmatrix} = V^{-1} \begin{pmatrix} X_0 \\ X_1 \\ \vdots \\ X_{p-1} \end{pmatrix} $$
where $V$ is the same matrix as in section 2 in this article, and $g$ is a linear function of $e_0,\ldots, e_n$ with some combination of the coefficients $a_1, \ldots, a_p$. For instance, if $p = 1$, we have $g(e_0,\ldots,e_n) = \sum_{k=0}^{n-1}a_1^k\cdot e_{n-k}$.
For an arbitrary $p$, the $k$-th term in the above sum (for the function $g$) is a polynomial of degree $k$, with $p$ variables $a_1, \ldots, a_p$. This formula also allows you to identify the auto-correlation structure, whether or not the time series is stationary or not.
Note
The matrix $V$ is a $p \times p$ Vandermonde matrix (its elements are the powers of the roots $r_1,\ldots,r_p$ of the characteristic polynomial $x^p = a_1 x^{p-1} + \cdots + a_{p-1} x + a_p$. The formula assumes that these roots are different (no multiple root) otherwise the formula must be adjusted.