function [M11,M2,FDD,F3,rc]=order2guts(BD,BDD,F1,F2,omega,sig);

% order2guts.m: CORRECTED CODE (APRIL 12, 2001)
% see order2test.m: tests of the code

%  function [M11,M2,FDD,F3,rc]=order2guts(BD,BDD,F1,F2,omega)     
%  System originally in the form Ka(z(t),z(t-1),eps(t),eta(t))=0, t=0,...,\infty, 
%  where Et[eta(t+1)]=0 for t\ge 0. Its second order expansion is K1 dz(t)= K2 
%  dz(t-1) + K3\eps(t) + K4\eta(t) + .5*(K11*dz(t)\otimes dz(t) + K12*dz(t)\otimes 
%  dz(t-1) ... It is assumed that the \eta terms enter linearly and thus generate no 
%  second-order terms. gensys transformations convert this form  to an equation 
%  where the K's are replaced by B's and z is are replaced by w=[y,x], with the 
%  following special characteristics:  the xy block of BD{1} is zero; the yy block of 
%  BD{1} is the identity; the y row of BD{4} is zero (and hence BD{4} is unused); the 
%  xx block of BD{2} is 
%  non-singular, and BD{2}{usix,usix}\BD{1}{usix,usix} (usix indexes the xx block of BD{2}) has 
%  all its eigenvalues 
%  less than (1/div)<1 in absolute value; the yy block of BD{1} has all its eigenvalues 
%  less than div>1 in absolute value. F1 and F2 are the coefficients in the first order
%  solution y(t)=F1*y(t-1)+F2*eps(t).  omega is the covariance matrix of eps.
%  Note that BD and BDD are cell arrays, whose elements are matrices.  BDD{j,k} for k<j
%  are not used, so should be left null to save space.

[nstate,nstate2]=size(F1); if nstate~=nstate2, error('F1 matrix not square'), end; 
[n,n2]=size(BD{1}); if n~=n2, error('BD{1} matrix not square'),end;   
uindex=nstate+1:n; 
sindex=1:nstate;
nu=n-nstate;
nerr=size(F2,2);
L=BD{2}(uindex,uindex)\BD{1}(uindex,uindex); 
M11 = -prodt(BDD{1,2}(uindex,sindex,sindex),F1,2,1); 

M11=M11+permute(M11,[1 3 2])-prodt(prodt(BDD{1,1}(uindex,sindex,sindex),F1,2,1),F1,2,1)...
   -BDD{2,2}(uindex,sindex,sindex);

% If do not use permute(M11,[1 3 2]): tHIS yields the quadratic form y'M11*y
% M11=2*M11-prodt(prodt(BDD{1,1}(uindex,sindex,sindex),F1,2,1),F1,2,1)...
%   -BDD{2,2}(uindex,sindex,sindex);

M11(:)=BD{2}(uindex,uindex)\M11(:,:);
%------------------Doubling algorithm-------------------------
R=F1; 
eest=1; 
Minc=zeros(size(M11)); 
crit= 10*sqrt((n-nstate)*nstate*nstate)*eps; 
dblinc=0;
while eest>crit & dblinc<200
   Minc(:)=L*M11(:,:);
   M11=M11+prodt(prodt(Minc,R,2,1),R,2,1);
   L=L*L;
   R=R*R;
   eest=sqrt(sum(sum(sum(Minc.*Minc))));
   dblinc=dblinc+1;
end
%fprintf(1,'dblinc=%d, eest=%g\n',dblinc,eest);
if dblinc>200
   fprintf(1,'Inaccuracy %g+%gi in order2guts.m\n ', [real(eest) imag(eest)])
   rc=1;
else
   rc=0;
end

% ==========================================================

worka=prodt(prodt(M11,F1,2,1),F1,2,1); 
   S1=worka;
worka=reshape(-BD{1}(sindex,uindex)*worka(:,:),[nstate,nstate,nstate]);
   S2=worka;   
workb=prodt(BDD{1,2}(sindex,sindex,sindex),F1,2,1);
workb=workb+permute(workb,[1 3 2]);
% ORIGINAL CODE
%FDD{1,1}=prodt(prodt(BDD{1,1}(sindex,sindex,sindex),F1,2,1),F1,2,1)+worka+workb...
%   +BDD{2,2}(sindex,sindex,sindex);

% CORRECTION, APRIL 12, 2001 
workc=prodt(BD{2}(sindex,uindex),M11,2,1);
FDD{1,1}=prodt(prodt(BDD{1,1}(sindex,sindex,sindex),F1,2,1),F1,2,1)+worka+workb...
   +workc + BDD{2,2}(sindex,sindex,sindex);



%----------------
worka=prodt(prodt(M11,F2,2,1,3),F2,2,1,3);
worka(:)=BD{1}(uindex,uindex)*worka(:,:);
workb=prodt(prodt(BDD{1,1}(uindex,sindex,sindex),F2,2,1,3),F2,2,1,3);
workc=prodt(BDD{1,3}(uindex,sindex,:),F2,2,1,3);
workc=workc+permute(workc,[1 3 2]);
worka=worka-workb-workc-BDD{3,3}(uindex,:,:);
M2=.5*((BD{2}(uindex,uindex)-BD{1}(uindex,uindex))\(worka(:,:)*omega(:)));
%----------------
worka=prodt(prodt(M11,F1,2,1),F2,2,1,3);
worka=reshape(-BD{1}(sindex,uindex)*worka(:,:),[nstate,nstate,nerr]);
  HK1=worka;
  
% THIS IS ORIGINAL CODE  
% workb=prodt(BDD{1,2}(sindex,sindex,sindex),F2,3 [!],1)...
% +permute(prodt(BDD{1,3}(sindex,sindex,:),F1,2,1,3),[1 3 2])+BDD{2,3}(sindex,sindex,:);

% APRIL 12, 2001: changed code: prodt(BDD{1,2}(sindex,sindex,sindex),F2, 2[!] ,1)...
% November 29, 2001 (CAS): To handle one-dimensional state vector
%  new code: prodt(BDD{1,2}(sindex,sindex,sindex),F2,2,1,3[!])...
workb=prodt(BDD{1,2}(sindex,sindex,sindex),F2,2,1,3)...
   +permute(prodt(BDD{1,3}(sindex,sindex,:),F1,2,1,3),[1 3 2])+BDD{2,3}(sindex,sindex,:);


HK2=prodt(prodt(BDD{1,1}(sindex,sindex,sindex),F1,2,1),F2,2,1,3);
HK3=prodt(BDD{1,2}(sindex,sindex,sindex),F2,2,1);
HK4=permute(prodt(BDD{1,3}(sindex,sindex,:),F1,2,1,3),[1 3 2]);
HK5=BDD{2,3}(sindex,sindex,:);

FDD{1,2}=prodt(prodt(BDD{1,1}(sindex,sindex,sindex),F1,2,1),F2,2,1,3)+worka+workb;


%--------------
worka=prodt(prodt(M11,F2,2,1,3),F2,2,1,3);
worka=reshape(-BD{1}(sindex,uindex)*worka(:,:),[nstate,nerr,nerr]);
worka=worka+prodt(prodt(BDD{1,1}(sindex,sindex,sindex),F2,2,1,3),F2,2,1,3);
workb=prodt(BDD{1,3}(sindex,sindex,:),F2,2,1,3);
workb=workb+permute(workb,[1 3 2]);
FDD{2,2}=worka+workb+BDD{3,3}(sindex,:,:);
%---------------------


% ORIGINAL CODE
% F3=-BD{1}(sindex,uindex)*M2;

% CORRECTION, APRIL 12, 2001
F3=(-BD{1}(sindex,uindex)+BD{2}(sindex,uindex))*M2; 


yxc=1001;
save order;

return;

% =============================================================================
% ====================== TESTS (see order2test.m) =============================
% 
% APRIL 12, 2001: TESTING ACCURACY OF M11 (test whether equation (16) holds)
LHS=prodt(prodt(prodt(BD{1,1}(uindex,uindex),M11,2,1,2),F1,2,1,3),F1,2,1,3);
RHS1=prodt(BD{1,2}(uindex,uindex),M11,2,1,2);
RHS2=prodt(prodt(BDD{1,1}(uindex,sindex,sindex),F1,2,1,3),F1,2,1,3);
RHS3=2*prodt(BDD{1,2}(uindex,sindex,sindex),F1,2,1,3);

RHS3=prodt(BDD{1,2}(uindex,sindex,sindex),F1,2,1,3);
RHS3=RHS3+permute(RHS3,[1,3,2]);
RHS4=BDD{2,2}(uindex,sindex,sindex);
test=-LHS+RHS1+RHS2+RHS3+RHS4;
if maxa(maxa(test))>1e-8; error('order2guts: computation of M11 inaccurate'); 
   end;
if maxaa(M11-permute(M11,[1 3 2])) > 1e-10; error('M11 not symmetric'); end;

% =============================================================================

G1=-BD{1}(sindex,uindex);   G2=BD{2}(sindex,:);       G3=BD{3}(sindex,:);
G11=BDD{1,1}(sindex,:,:);   G12=BDD{1,2}(sindex,:,:); G13=BDD{1,3}(sindex,:,:); 
                            G22=BDD{2,2}(sindex,:,:); G23=BDD{2,3}(sindex,:,:); 
                                                      G33=BDD{3,3}(sindex,:,:);
                                             
                                                            
% ====== USE SIMULATION TO TEST WHETHER EQUATIONS (10) and (12) are met
yL=randn(cols(sindex),1)*100;
yL=zeros(cols(sindex),1);
eH=randn(cols(F2),1)*100;
yH=F1*yL+F2*eH;       
x=zeros(cols(uindex),1);

sig=1; % value of 'sig' set here is irrelevant (provided sig>0)

R1=G1*(.5*prodt(prodt(M11,yH,2,1),yH,2,1)+M2*sig^2);  

xL=.5*prodt(prodt(M11,yL,2,1),yL,2,1)+M2*sig^2;
R2=G2*[yL;xL];
R3=G3*eH;

R4=.5*prodt(prodt(G11,[yH;x],2,1),[yH;x],2,1);
R5=prodt(prodt(G12,[yH;x],2,1),[yL;x],2,1);
R6=prodt(prodt(G13,[yH;x],2,1),eH,2,1);
R7=.5*prodt(prodt(G22,[yL;x],2,1),[yL;x],2,1);
R8=prodt(prodt(G23,[yL;x],2,1),eH,2,1);
R9=.5*prodt(prodt(G33,eH,2,1),eH,2,1);
TY1=R1+R2+R3+R4+R5+R6+R7+R8+R9;  %value of dy(t) determined by eqn (10)

% &&&&&&&&&&&&&
% NEXT LOOK AT EQUATION (12)

S1=F1*yL+F2*eH+F3*sig^2;
S2=.5*prodt(prodt(FDD{1,1},yL,2,1),yL,2,1);
S3=prodt(prodt(FDD{1,2},yL,2,1),eH,2,1);
S4=.5*prodt(prodt(FDD{2,2},eH,2,1),eH,2,1);
   
TY2=S1+S2+S3+S4;  % value of dy(t) determined by eqn (12)
test=TY1-TY2;                                                                                                                                                                             

if maxa(test) > 1e-8; error('order2guts: inaccurate solution'); end;

save order;