An AR($p$) process is any causal and weakly stationary solution to the equations $$ X_t = \beta_1 X_{t-1} + \dotsc + \beta_p X_{t-p} + \epsilon_t, \qquad t \in \mathbb{Z} $$ where the polynomial $B(z) = 1 - \beta_1z - \dotsc - \beta_p z^p$ must not have roots on the unit disk. Here, I assume the noise $\epsilon_t$ to be iid with $E\epsilon_t = 0, \text{Var}(\epsilon_t) = \sigma^2$.
There is a unique solution to these equations which fulfills weak stationarity and causality. It can be written as a MA($\infty$): $$ X_t = \sum_{k = 0}^\infty \chi_k \epsilon_{t-k} $$ with appropiate coefficients $\chi_k$.
I believe that $X_t$ is a strictly stationary process due to the $\epsilon_t$ being iid, but in the literature I can't find a proof for this claim (all proofs that I saw are about weak stationarity).
How does the proof work?
EDIT: Note that I am not fixing a "starting distribution" for $X_0$. This is unnecessary because the process starts at time $t = -\infty$.