Generating autocorrelated random numbers with a specific distribution I want to generate autocorrelated data that follows a standard normal distribution. Following the algorithm in this paper :
$$ X_t= \frac{\alpha X_{t-1}+(1-\alpha)Z_t}{\sqrt{(1-\alpha)^2+\alpha^2}}$$
$$X_1=Z_1$$
where $Z_t$ is normally distributed with 0 mean and 1 variance. The variance of $X_t$ is given by
$$Var(X_t)=1 \text{   for   } t=1 $$
$$ Var(X_t)=\frac{(1-\alpha)^2Var(Z_t)}{(1-\alpha)^2+\alpha^2}+\frac{\alpha^2Var(X_{t-1})}{(1-\alpha)^2+\alpha^2} \text{   for   } t>1 $$
Because $Var(X_1)$ is 1, $Var (X_2)$ can be shown to be 1 and thus $Var(X_3)=1$ and so on,hence
$$X_t~N(0,1)\forall t$$
I understand everything so far. However, when I implement this recursive algorithm in Matlab my $X_t$ does not have 0 mean and 1 standard deviation even though I generate the $Z_t$ values from a standard normal distribution. The distortion becomes huge at some $\alpha$ values. Can someone shed some light on why my $X_t$ does not follow a standard normal distribution? Am I missing a condition in the above equations?
Thanks

The sequence $(X_t)$ is a Markov chain
$$X_t=\beta X_{t-1}+\sigma\epsilon_t\qquad\beta^2+\sigma^2=1$$
and the terms are correlated with
$$\text{cor}(X_t,X_{t+k})=\beta^k$$
This correlation implies that the convergence of the empirical moments is slower than with an iid sequence.
For instance here is a basic R code reproducing the generation of the sequence
simar=function(al=.5,omal=1-al,N=1e2){
x=z=rnorm(N)
for(t in 2:N)x[t]=(al*x[t-1]+omal*z[t])/sqrt(al^2+omal^2)
print(c(mean(x),var(x)))}

Then
> simar()
[1] -0.2406807  0.9394529
> simar()
[1] 0.05408356 1.55802573
> simar(N=1e4)
[1] 0.006863719 1.007094258
> simar(N=1e6)
[1] 0.001068688 1.003048360

reflects on the slow convergence. The closer $\beta$ is to $1$, the slower the convergence, e.g.
> simar(al=.9)
[1] 1.687833 0.187142
> simar(al=.9)
[1] 1.32906836 0.07415279
> simar(al=.9,N=1e6)
[1] -0.01555348  0.95794239
> simar(al=.9,N=1e6)
[1] -0.002255304  1.000054363

