Is the AR(1) process always Gaussian given Gaussian innovations? I found that $AR(1)$ process $x_t=\phi x_{t-1}+\epsilon_t$ is not always Gaussian given Gaussian innovations $\epsilon_t$. This only happens when the $AR(1)$ model coefficient is very large. This goes against with the theory (If the white noise $\epsilon_t$ is a Gaussian process then $x_t$ is also a Gaussian process). How to explain this then? The following R codes gives that only $363$ out of the 1000 generations of $AR(1$) processes are gaussian process by the  Shapiro-Wilk  test.
decision=c()
for (k in 1:1000){ # 1000 runs 
  phi=0.9 # AR(1) coefficient
  x=0.5 # value of x_0
  X=c() # generated AR(1) process
  for (i in 1:1000){  # length of the process is 1000
    x=phi*x+rnorm(1) # Gaussian innovations
    X=c(X,x) 
  }  
  decision=c(decision,shapiro.test(X)$p.value) # Gaussian test
}
sum(decision>0.05)/1000 # percentage of Gaussian processes

 A: By way of summary of the comments, so that this question has an answer (if Alecos would like to present a summary I'll happily delete this and upvote it instead, as long as the question ends up with an answer) --
As Alecos points out, the AR(1) process for $\phi=0.9$ is definitely Gaussian.
Since it's stationary, observations from it should all have the same Gaussian distribution; the OP is right to expect that it should, since it's undeniably so. As Alecos says, the pen-and-paper result should suggest there's a problem with either the code or the test.
There are in fact several such problems:
(1) The Shapiro Wilk test assumes independence. We don't have it so the test doesn't apply. The OP notes that the test doesn't seem to have problems even with moderately large values of the parameter, and that's not surprising - because of the way the test works, it should be fairly robust to mild dependence. The Shapiro-Wilk doesn't look at the dependence in consecutive values. The most noticeable effect of the dependence on the distribution of the order statistics will be to increase their variance, but the Shapiro Wilk won't notice that at all. There will be a tendency for the tails to wander more from the straight line than with an independent series and eventually that sort of deviation will become detectable.
(2) The code doesn't have any warmup period. The mean and variance of the early values won't correspond to the mean and variance of the AR(1) process (so we don't have independence OR identically distributed values).
