%This function implements the method presented in Tauchen (1986) to approximate a continuous AR(1) process of the form: zt+1=rho zt +epsilont+1.
%The code requires as inputs the standard deviation for the error term (sigmae),
%persistence parameter (rho), number of states (m), and grid size (lambda).
function[P,z]=tauchenAR(rho,sigmae,lambda,m)
z1=-(lambda*sigmae)/(sqrt(1-rho^2));
step=-2*z1/(m-1);
for i=1:m
    z(i,1)=z1+(i-1)*step;
end
for i=1:m
    pi1(i,1)=normcdf(((z1-rho*z(i))/sigmae)+(step/(2*sigmae)));
    for j=2:m-1
        pij(i,j)=normcdf(((z(j)-rho*z(i))/sigmae)+(step/(2*sigmae)))-normcdf(((z(j)-rho*z(i))/sigmae)-(step/(2*sigmae)));
    end
end
pij(:,1)=pi1;
pim=1-sum(pij,2);
P=[pij pim];

        