function [y,P,pi] = tauchen(lambda,sige,N,m,nquad)
%TAUCHEN   finds the transition probabilities for a Markov chain intended
%          to mimic the properties of the autoregressive process:
%
%               y[t] = lambda y[t-1] + e[t],  e[t] ~ N(0,sige^2).     (1)
%
%          [y,P,pi]=tauchen(lambda,sige,N,m,nquad) constructs the transition 
%          matrix P and stationary distribution pi from equation (1) for 
%          an N-state Markov-chain, y=[y1,..yN].  The value yN is a multiple 
%          m of the unconditional standard deviation for y[t], e.g., 
%          sigy=sige/sqrt(1-lambda^2).  nquad is the number of quadrature
%          points when integrating over the normal c.d.f.  If only three 
%          inputs are provided, i.e., tauchen(lambda,sige,N), then m is 
%          set equal to 3 and nquad is set equal to 20.

%          References: 
%           [1] Tauchen, G. ``Finite State Markov-chain approximations
%               to Univariate and Vector Autoregressions,'' Economics 
%               Letters 20, pp. 177-181, 1986.

%          Ellen McGrattan, 9-19-93
%          Revised, ERM, 11-2-93

if nargin<4; 
  m     = 3; 
  nquad = 20; 
end;
sigy    = sige/sqrt(1-lambda*lambda);
y(N)    = m*sigy;
y(1)    = -y(N);
w       = 2*y(N)/(N-1);
y(2:N-1)= y(1)+w:w:y(N)-w;
minf    = (y(1)-lambda*y(N)+0.5*w)/sige-sigy;   
for j=1:N;
  a              = minf;
  b              = (y(1)-lambda*y(j)+0.5*w)/sige; 
  p              = zeros(nquad+1);
  p(1,nquad+1)   = 1; 
  p(2,nquad)     = 1;
  for l=2:nquad;  
    p(l+1,:)     = ( (2*l-1)*[p(l,2:nquad+1),0]-(l-1)*p(l-1,1:nquad+1) )/l;
  end;
  x              = sort(roots(p(nquad+1,:)));
  A              = p(1:nquad,2:nquad+1)*( (ones(nquad,1)*x').^ ...
                   ([nquad-1:-1:0]'*ones(1,nquad)) );
  wgt            = A\[b-a;zeros(nquad-1,1)];
  x              = sort(.5*(a*(1+x)+b*(1-x)));
  P(j,1)         = 3.989422804014327e-01*wgt'*exp(-x.*x/2);
  for k=2:N-1;
    a            = (y(k)-lambda*y(j)-0.5*w)/sige;
    b            = (y(k)-lambda*y(j)+0.5*w)/sige;
    p            = zeros(nquad+1);
    p(1,nquad+1) = 1; 
    p(2,nquad)   = 1;
    for l=2:nquad;  
      p(l+1,:)   = ( (2*l-1)*[p(l,2:nquad+1),0]-(l-1)*p(l-1,1:nquad+1) )/l;
    end;
    x            = sort(roots(p(nquad+1,:)));
    A            = p(1:nquad,2:nquad+1)*( (ones(nquad,1)*x').^ ...
                  ([nquad-1:-1:0]'*ones(1,nquad)) );
    wgt          = A\[b-a;zeros(nquad-1,1)];
    x            = sort(.5*(a*(1+x)+b*(1-x)));
    P(j,k)       = 3.989422804014327e-01*wgt'*exp(-x.*x/2);
  end;
  P(j,N)         = 1-sum(P(j,1:N-1));
end;
if nargout>1;
  pi    = ones(N,1)/N;
  for i = 1:100;
    pi  = (pi'*P)';
  end;
end;


