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;