% % % We seek a solution of the type
% % 
% % c=Pi21*s+PHI2*x
% % 
% % s=Pi11*s(-1)+PHI1*x(-1);
% 
% 
namemod=input('give name of model: ','s');

eval(['load matrices_',namemod,]);
% 
% %%%%%%%%%%%%%%%%%%%%%%%%5
% nolev=matrices.nolev;
% jumping=matrices.jumping;
% 
% shocky=matrices.shocky;
% varyt=matrices.varyt;
% nsh=size(shocky,1);
% prede=matrices.neq-jumping;
% neq=size(varyt,1);
% 
% mxpl=matrices.mxpl;
% mxmi=matrices.mxmi;
% 
% M_=struct('endo_names',varyt(end:-1:1,:));
% M_.fname=namemod;
% M_.param_names=matrices.parnames;
% 
% lgy_=M_.endo_names;
% assignin('base','lgy',lgy_);
% assignin('base' ,'M_',M_);
% % keyboard
% if exist([namemod,'_steadystate.m'],'file')
%     ys=feval([namemod,'_steadystate']);
%     for jj=1:length(ys);
%         eval([lgy_(jj,:),'=',num2str(ys(jj),16),';']);
%                 eval(['global  ',deblank(lgy_(jj,:)),'_ss']);
% % Assume that the ss value of the variables is computed in the steadystate file and the _ss are assigned there
% %         eval([deblank(lgy_(jj,:)),'_ss=',num2str(ys(jj),16),';']);
%     end
%     M_=evalin('base','M_');
% %      keyboard
%     NMS=str2cell(M_.param_names);
%     for jj=1:length(M_.params);
%         eval([M_.param_names((jj),:),'=',num2str(M_.params(jj),16),';']);
%     end
% else % it is assumed that the parameter values are assigned in the steady_state file
%     eval(['parameters_',namemod]);
% tmp=[];
% for jj=1:size(matrices.parnames,1);
%     tmp(jj,1)=eval([matrices.parnames(jj,:)]);
% end
% M_.params=tmp;
% end
% %%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% param_neoc_grwt_dynare;
% 
% % Get the canonical form
% F1=eval(matrices.A_plus1);
% F2=eval(matrices.A_0);
% F3=eval(matrices.A_minus1);
% % model is F1*x(+1)+F2*x+F3*x(-1)=0; need to change sign of matrices
% % accordingly
% 
% 
% %%% eliminate the shock from this form: these take the last nsh
% %%% rows/columns
% tmp=F1(1:end-nsh,end-nsh+1:end);
% F1=F1(1:end-nsh,1:end-nsh);
% eval([namemod,'_X']);
% N=Xmx;
% % N=F3(end-nsh+1:end,end-nsh+1:end);
% tmp2=F3(1:end-nsh,end-nsh+1:end);
% tmp3=sum(abs(tmp2),2);% find which shocks enter the model with t-1
% tmp4=F2(1:end-nsh,end-nsh+1:end);
% F2=F2(1:end-nsh,1:end-nsh);
% F3=F3(1:end-nsh,1:end-nsh);
% [jj]=find(tmp3>0);
% if ~isempty(jj);
%     F1=[F1,zeros(size(F1,1),length(jj));
%         zeros(length(jj),size(F1,2)+length(jj))];
%     F2=[F2,tmp2(:,jj);% add the pointer to lagged shocks as pointer to new variables
%         zeros(length(jj),size(F1,2)),eye(length(jj))];
%     F3=[F3,zeros(size(F1,1),length(jj));
%         zeros(length(jj),size(F1,2)+length(jj))];                
% 
% end
%        C=zeros(size(F3,1),nsh);
%     C(end-length(jj)+1:end,jj)=-1; 
% C(1:size(tmp4,1),:)=C(1:size(tmp4,1),:)-tmp*N+tmp4;% point to time t shock;
% 
% 
% 
% neq=size(matrices.A_0,1)-nsh+length(jj);
% 
% jumping=size(F1);
% prede=neq+length(jj);
% 
% 
% 
% % 
% % % write the system in companion form
% % 
% 
% AA=[-F1,zeros(size(F1));
%     zeros(size(F1)),eye(size(F1,1))];
% BB=[F2,F3;
%     eye(size(F1,1)),zeros(size(F2))];
% 
% C=[C;
%     zeros(size(F1,1),size(C,2))];
% 
% 
% % %AA=eval(matrices.A);
% % eval([namemod,'_A']);
% % % take out the shock part (end of matrix)
% % AA=AA(1:end-nsh,1:end-nsh);
% % 
% % %BB=eval(matrices.B);
% % eval([namemod,'_B']);
% % % take out shock part
% % C=BB(1:end-nsh,end-nsh+1:end);
% % BB=BB(1:end-nsh,1:end-nsh);
% % %N=eval(matrices.Xmx);
% % eval([namemod,'_X']);
% % N=Xmx;
% % 
% % neq=neq-nsh;
% % jumping=jumping-nsh;
% % C=zeros(size(AA,1),nsh);%set dummy shocks to 2 as 1 could generate problems
% % C((mxpl-1)*neq+neq-nsh+1:(mxpl-1)*neq+neq,1:nsh)=eye(nsh);%
% 
% AA=AA(:,[end:-1:1]); % jump last
% BB=BB(:,[end:-1:1]);
% 
% [s,t,q,z]=qzf(AA,BB,1);
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%5  MUC
% 
% s11=s(1:prede,1:prede);
% s12=s(1:prede,prede+1:end);
% s22=s(prede+1:end,prede+1:end);
% 
% t11=t(1:prede,1:prede);
% t12=t(1:prede,prede+1:end);
% t22=t(prede+1:end,prede+1:end);
% 
% zt=z';% its' already transposed
% 
% 
% zt11=zt(1:prede,1:prede);
% zt12=zt(1:prede,prede+1:end);
% zt21=zt(prede+1:end,1:prede);
% zt22=zt(prede+1:end,prede+1:end);
% 
% z11=z(1:prede,1:prede);
% z12=z(1:prede,prede+1:end);
% z21=z(prede+1:end,1:prede);
% z22=z(prede+1:end,prede+1:end);
% 
% Chat=q'*C;
% 
% %%%%%%%%%%%%%%%%
% % CHECK OUT THE SIGN
% Pi21=-zt22\zt21;
% Pi11=z11*(s11\t11)/z11;
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%
% 
% G1=(s11*(zt11+zt12*Pi21));
% G2=s11*zt12+s12*zt22;
% G3=t11*zt12+t12*zt22;
% 
% PHI2=zeros(prede,nsh);
% if size(s22,1)*nsh<200;
%     
%     tmp=-(kron_d(N',s22*zt22)+kron_d(eye(nsh),t22*zt22))\Chat(prede+1:end,:);
%     
%     
%     
%     PHI2=reshape(tmp,prede,nsh);
% else
%     tmp=inv(t22*zt22);
%     %
%     PHI2_=zeros(size(PHI2));
%     for jj=1:500;
%         PHI2_=-tmp*(Chat(prede+1:end,:)+(s22*zt22*PHI2_*N));
%     end
% end
% 
% PHI1=G1\(G2*PHI2*N+G3*PHI2+Chat(1:prede,:));

a=cputime;
[Pi11,Pi21,PHI1,PHI2]=solver_muc(namemod,matrices);
d=cputime;
[Pi,PHI,AA_,BB_,numpre,aqz,bqz,Qqz,Zqz,inz22,algo]=firstorder_stuff(namemod,matrices,0);
e=cputime;
sprintf('MUC %2.5f Klein %2.5f',d-a,e-d)
%%%%%% test result

shock=1;
tmp1=0*PHI1*shock;
tmp2=[shock;zeros(3,1)];
for kk=2:10
    shock(1,kk)=N*shock(1,kk-1);
    tmp1(:,kk)=Pi11*tmp1(:,kk-1)+PHI1*shock(1,kk-1);
    tmp2(:,kk)=PHI*tmp2(:,kk-1);
end
disp('discrepancy')
[Pi(2:end,:)*tmp2-(Pi21*tmp1+PHI2*shock)]