Weakly stationary Gaussian AR(1) process is strict stationary? 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.
 A: There is more to a Gaussian process than just all the random variables in the process being Gaussian random variables individually: it is necessary that every (finite) collection of (two or more) random variables from the process also be jointly Gaussian or have a multivariate normal distribution depending on which flavor of wording you prefer.  For a large collection of ways in which individually Gaussian random variables can fail to be jointly Gaussian, see here.
However, in this instance, what you are saying is correct. A weakly stationary Gaussian process is also a strictly stationary process because a multivariate normal distribution depends only on the means (known to be constant), variances (also known to be constant), and covariances of the random variables involved, and the covariances depend only on the time differences between the variables. Thus,
$\operatorname{cov}(x_i, x_j) = \operatorname{cov}(x_{i+m}, x_{j+m})$ and so the multivariate normal distribution of $(x_i, x_j, x_k, \cdots)$ is the same as the multivariate normal distribution of $(x_{i+m}, x_{j+m}, x_{k+m}, \cdots)$ which is what we need in order to assert strict stationarity.
A: 
I know the process is weakly stationary.

No, you don't.  In fact, just as with establishing strong stationarity, in order to get weak stationarity you need to make an additional assumption about an "anchoring point".
It is a common misconception in time-series analysis to think that a root inside the unit circle is a sufficient condition for stationarity (see here for further discussion).  The recursive equation is insufficient to determine stationarity without setting a "starting distribution" for an anchoring point for the process.  To see that what you have is insufficient, consider the following $\text{AR}(1)$ process that satisfies your recursive equation and error distribution, but is not even weakly stationary.

An example with $|a|<1$ that is not even weakly stationary: Suppose you consider a time-series model that satisfies your $\text{AR}(1)$ recursive equation and error distribution, and you specify a "starting distribution" for the anchoring point $X_0$ that is still a normal distribution, but with some arbitrary mean and variance:
$$X_0 \sim \text{N} ( \mu_0, \sigma_0^2 ).$$
(Although I omit it from the notation, strictly speaking this is the conditional distribution, conditioning on the entire error series ---i.e., the anchoring point is assumed to be independent of the error series.)  With this "starting distribution" it can be shown that the series of marginal distributions is:
$$X_t \sim \text{N} \Big( a^t \mu_0, a^{2t} \sigma_0^2 +\frac{1-a^{2t}}{1-a^2} \sigma_e^2 \Big)
\quad \quad \quad \text{for } t \in \mathbb{Z}.$$
For $\mu_0 \neq 0$ this series has non-stationary mean, though it converges to zero mean.  For $\sigma_0 \neq \sigma_e$ this series has non-stationary variance, though it converges to the asymptotic variance $\sigma_e^2 / (1-\alpha^2)$.  In this case the process is not even first-order stationary.
The simplest example of this is to take $\sigma_e = \sigma_0 = 0$ so that the anchoring point is the fixed point $\mu_0 \neq 0$ and the error terms in the series are all zero.  In this case we have the non-statinary deterministic sequence $X_t = a^t \mu_0$ for all $t \in \mathbb{Z}$.  This sequence satisfies recursive equation $(1)$ but it has a non-constant mean.

Establishing weak and strong stationarity: For the $\text{AR}(1)$ process defined by your recursive equation and error distribution you need two things to establish stationarity.  You need to have $|a|<1$ and you need to have a "starting distribution" for an anchoring point that ensures the required stationarity condition (and which implicitly assumes that the anchoring point is independent of the error series).  To establish weak stationarity in your process (i.e., stationarity of mean and variance) you need to have a starting distribution with mean $\mu_0=0$ and variance $\sigma_0^2 = \sigma_e^2 / (1-\alpha^2)$.  To establish strict stationarity you need the additional assumption that this starting distribution is Gaussian.
To extend this requirement to broader $\text{ARMA}$ models, you need to impose $p$ anchoring points, where you have $p$ auto-regressive terms in your model.  Thus, for example, establishing stationarity in an $\text{ARMA}(2,3)$ model would require a starting distribution for two anchoring points.
