rho = input('What value for rho? (.95^4, .85^4, .75^4, 0) ');
if rho==.95^4;
  sig = .0136;
elseif rho==.85^4;
  sig = .02;
elseif rho==.75^4;
  sig = .02223;
else
  sig = .02345;
end;

lz1 = log(2.16); 
lz2 = log(1.54); 
nz  = 11;

[y,P,pi]=tauchen(rho,sig,nz);
del = y(2)-y(1);
z1  = exp(lz1+[y([6,4,3,4,6,8,10,11]),y(11)+.5*del,y(11)+.5*del, ...
               y([11,10,8,6,6]),y]);
[y,P,pi]=tauchen(rho,sig*2,nz);
del = y(2)-y(1);
z2  = exp(lz2+[y([6,4,3,4,6,8,10,11]),y(11)+.5*del,y(11)+.5*del, ...
               y([11,10,8,6,6]),y]);

Pp  = P';
Pp  = [Pp(:,[6,4,2,4,6,8,10,11,11,11,11,10,8,6,6]),Pp];

disp('Labor-Augmenting Technologies:')
disp([z1;z2])

disp('Transition matrix:')
disp([zeros(15,26);Pp])
%disp([diag(ones(14,1),-1)+diag([zeros(14,1);1]),zeros(15,11); zeros(11,15),P'])


