Let the sequence of random variables $x_t$ be given by $$x_t=ax_{t-1}+e_t,\quad t\in\mathbb{Z},\quad (1)$$ where $e_t\sim i.i.d.N(0,\sigma_e^2)$, $\lvert a\rvert\leq 1$. I know the process is weakly stationary. So $E(x_t)=E(x_{t-1})=0$ and $Var(x_t)=Var(x_{t-1})=\sigma_e^2/(1-a^2)$ for any $t$. Also, $x_{t-1}=\sum_{j=0}^\infty a^je_{t-1-j}$ is gaussian and independent from $e_t$, which is also gaussian. Then $x_t$ must be gaussian.
These observations, seems to confirm that the $AR(1)$ process is strictly stationary. If I define (1) in terms of times series for the indices $t\in\{1,2,3,\dotsc\}$, perhaps some additional assumption is needed, e.g., $x_0\lvert (e_1,e_2,\dotsc) \sim N(0,\sigma_e^2/(1-a^2))$.
Question: Are these observations correct? Do you see any flaw?
Comments
I will sketch my arguments against the answer below.
In Ben's answer, it was said that model (1) does not imply weak stationarity. His arguments are based on model (1) defined on $t\geq0$. In model (1), I assumed $t\in\mathbb{Z}$ and for any integer $t$, the error term is gaussian. Writing $x_{t}=\sum_{j=0}^\infty a^j e_{t-1-j}$ we see that $E(x_t)=0$ for any integer $t$, and since $\sum_{j=0}^\infty a^{2j}\sigma_e^2=\sigma_e^2/(1-a^2)<\infty$, $Var(x_t)=\sigma_e^2/(1-a^2)$ for any integer $t$, as well.
Th
For the covariance, \begin{align} Cov(x_t,x_{t+k})&=Cov(\sum_ja^j e_{t-j}, \sum_l a^l e_{t+k-l})=\sum_{j,l}a^{j+l}Cov(e_{t-j},e_{t+k-l})\\ &=\sum_{j=0}^\infty \sum_{w=-k}^\infty a^{j+w+k} Cov(e_{t-j},e_{t-w})=\sigma_e^2\sum_{j=0}^\infty a^{2j+k}\\ &=\sigma_e^2a^k/(1-a^2). \end{align}
These suggests that I do know that the process is weakly stationary.