function [yd,P,pi,psum,S] = tauch2D(A,Sige,N,grid,yobs)
%TAUCH2D     finds the transition probabilities for a Markov chain intended
%            to mimic the properties of the autoregressive process for 
%            the 2x1 vector y:
%
%            (1)  y[t] = A y[t-1] + e[t],  e[t] ~ N(0,Sige), Sige not diagonal.
%
%            [y,P,pi,psum,S]=tauch2D(A,Sige,N,grid,yobs) computes 
%            the transition matrix P and stationary distribution pi from
%            equation (1) for an N*-state Markov-chain (N*=N(1)N(2)).
%            The ith value for the first variable in y is the midpoint
%            of [grid(1,i),grid(1,i+1)].  Similarly the ith value for the
%            second variable is the midpoint of [grid(2,i),grid(2,i+1)].
%            If the user has time series for y[t], namely the (nt x 2) 
%            matrix yobs, then the program computes the indices S for 
%            the Markov Chain that most closely mimic the observed 
%            series.  Output psum is a check on the domain sizes -- psum 
%            is the sum of rows in P before they are normalized to sum 
%            to 1.  If only three inputs are provided, then even spacing
%            is used for all variables with the first node at -3 times
%            the standard deviation of the variable and the last at +3
%            times the standard deviation.  With no yobs input, S=0.
%
%            See also VECTAUCH.m which assumes Sige is diagonal MxM
%            as in Tauchen's article.  See also ESTAUCH2D.m which
%            uses even spacing by default.

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

%            Ellen McGrattan, 1-21-97
%            Revised, ERM, 6-19-03

[M,tem]        = size(A);
if nargin<4; m = 3; end;
if M~=2;
  error('TAUCH2D for cases with 2-dimensional system')
end;
N              = N(:);
Sigy           = dlyap(A,Sige);
sigy           = sqrt(diag(Sigy));
TN             = prod(N);
y              = zeros(2,max(N));
wlag           = y;
wlead          = y;
for i=1:2;
  if nargin<4;
    y(i,1:N(i))     = linspace(-3*sigy(i),3*sigy(i),N(i));
    wlag(i,1:N(i))  = y(i,N(i))/(N(i)-1);
    wlead(i,1:N(i)) = wlag(i,1:N(i));
  else
    y(i,1:N(i))     = .5*grid(i,1:N(i))+.5*grid(i,2:N(i)+1);
    wlag(i,1:N(i))  = y(i,1:N(i))-grid(i,1:N(i));
    wlead(i,1:N(i)) = grid(i,2:N(i)+1)-y(i,1:N(i));
  end
end;

%
% ind = indices of the TN discrete states
%
yd             = zeros(TN,2);
wm             = zeros(TN,2);
wp             = zeros(TN,2);
ind            = ones(TN,2);
ind(:,1)       = vec([1:N(1)]'*ones(1,N(2)));
ind(:,2)       = vec(ones(N(1),1)*[1:N(2)]);
yd(:,1)        = y(1,ind(:,1))';
yd(:,2)        = y(2,ind(:,2))';
wm(:,1)        = wlag(1,ind(:,1))';
wm(:,2)        = wlag(2,ind(:,2))';
wp(:,1)        = wlead(1,ind(:,1))';
wp(:,2)        = wlead(2,ind(:,2))';

P              = zeros(TN);
for j=1:TN;
  mu           = A*yd(j,:)';
  for l=1:TN;
    a1         = yd(l,1)-mu(1)-wm(l,1);
    b1         = yd(l,1)-mu(1)+wp(l,1);
    a2         = yd(l,2)-mu(2)-wm(l,2);
    b2         = yd(l,2)-mu(2)+wp(l,2);
    [x1,wgt1]  = qgausl(a1,b1,20);
    [x2,wgt2]  = qgausl(a2,b2,20);
    sum1       = 0;
    for i1=1:20
      for i2=1:20
        x      = [x1(i1);x2(i2)];
        sum1   = sum1+wgt1(i1)*wgt2(i2)*exp(-.5*x'*inv(Sige)*x);
      end;
    end;
    P(j,l)     = sum1/(6.28318530717959*sqrt(det(Sige)));
  end;
end;
P0=P;
psum           = sum(P');
P              = P./(psum'*ones(1,TN));
[V,D]          = eig(P');
tem            = find(diag(D)==max(diag(D)));
pi             = V(:,tem)/sum(V(:,tem));     % double check: pi = (pi'*P)'?

S              = 0;
if nargin>4;
  nt           = length(yobs);
  SS           = zeros(nt,2);
  for i=1:2;
    bins       = [grid(i,1:N(i))',grid(i,2:N(i)+1)'];
    bins(1,1)  = -inf;
    bins(N(i),2)= inf;
    for j=1:nt;
      find( yobs(j,i)> bins(:,1) & yobs(j,i) <= bins(:,2));
      SS(j,i)  = find( yobs(j,i)> bins(:,1) & yobs(j,i) <= bins(:,2) );
    end;
  end;
  for i=1:nt;
    k = 1;
    for j=1:TN;
      if (ind(j,:)==SS(i,:)); k=j; end;
    end;
    S(i,1) = k;
  end;
end;


