Generating Auto-Correlated Random Numbers with finite correlation length How can I generate a lot of (for example N=20000 or more) auto-correlated random numbers {x_n, n=1,2,...,N} with uniform or normal probability distribution and finite correlation length (Here eq. correlation length = 3), i.e. 
=a0*d(n-m,0)+a1*d(n-m,1)+a2*d(n-m,2)+a3*d(n-m,3)
where 1>=a0>=a1>=a2>=a3>=0 will given as inputs.
I write a code according to "Method 1 of simulation", but it takes a lot of Ram and time. 
Can any one write an optimized code for me about this problem in fortran, matlab or mathematica or propose a revised or new algorithm for it?
 A: I find the answer to my question. The answer is based on Davies and Harte method and generate n=2^g random numbers. Here sigma(1)=a0, sigma(2)=a1, and .... The matlab code is :
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g=25;
n=2^g;
sigma=[1 0.5];
mv=size(sigma); m=mv(2);
r=zeros(1,2*n);
for i=1:m
r(i)=sigma(i);

if (i < m)
    r(2*n-i+1)=sigma(i+1);
end

end
%disp(r);
lambda=fft(r);
%disp(lambda);
z=randn(1,2*n);
%disp(z);
y=zeros(1,2*n);
y(1)=sqrt(2*n*lambda(1))*z(1);
y(n+1)=sqrt(2*n*lambda(n+1))*z(2*n);
for i=1:n-1
y(i+1)=sqrt(n*lambda(i+1))*(z(2*i)+1i*z(2*i+1));
y(i+1+n)=sqrt(n*lambda(i+1+n))*(z(2*n-2*i)-1i*z(2*n-2*i+1));

end
%disp(y);
x=ifft(y);
%disp(x);
fileID = fopen('CorrelatedRandomNumbers.txt','w');
fprintf(fileID,'%d\n',x);
sizx=size(x);
disp([' nt = ' num2str(sizx(2)/2)]);
[ACF, lags, bounds] = autocorr(x, [], 2); 
fileID = fopen('AutoCorrelation.txt','w');
fprintf(fileID,'%d\n',ACF);
disp('End of Program ...')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
