PI_gensys

 System given as
        a0*E_t[y(t+1])+a1*y(t)=a2*y(t-1)+c+psi*eps(t)
 with z an exogenous variable process.
 Returned system is
       [s(t)' x(t)' E_t x(t+1)']'=G1pi [s(t-1)' x(t-1)' x(t)]'+C+impact*eps(t),
  and (a) the matrix nmat satisfying   nmat*E_t z(t)+ E_t x(t+1)=0
      (b) matrices TT1, TT2  that relate y(t) to these states: y(t)=[TT1 TT2][s(t)' x(t)']'.
 Note that the dimension of the state vector = dim(a0)+NO_FL_EQS

 If div is omitted from argument list, a div>1 is calculated.
 eu(1)=1 for existence, eu(2)=1 for uniqueness.  eu(1)=-1 for
 existence only with not-s.c. z; eu=[-2,-2] for coincident zeros.
 Based on
 Christopher A. Sims
 Corrected 10/28/96 by CAS

0001 function [G1pi,C,impact,nmat,TT1,TT2,gev,eu, DD, E2, E5, GAMMA, FL_RANK ]=PI_gensys(a0,a1,a2,a3,c,PSI,NX,NETA,FL_RANK,M_,options_)
0002 % System given as
0003 %        a0*E_t[y(t+1])+a1*y(t)=a2*y(t-1)+c+psi*eps(t)
0004 % with z an exogenous variable process.
0005 % Returned system is
0006 %       [s(t)' x(t)' E_t x(t+1)']'=G1pi [s(t-1)' x(t-1)' x(t)]'+C+impact*eps(t),
0007 %  and (a) the matrix nmat satisfying   nmat*E_t z(t)+ E_t x(t+1)=0
0008 %      (b) matrices TT1, TT2  that relate y(t) to these states: y(t)=[TT1 TT2][s(t)' x(t)']'.
0009 % Note that the dimension of the state vector = dim(a0)+NO_FL_EQS
0010 %
0011 % If div is omitted from argument list, a div>1 is calculated.
0012 % eu(1)=1 for existence, eu(2)=1 for uniqueness.  eu(1)=-1 for
0013 % existence only with not-s.c. z; eu=[-2,-2] for coincident zeros.
0014 % Based on
0015 % Christopher A. Sims
0016 % Corrected 10/28/96 by CAS
0017 
0018 % Copyright (C) 1996-2009 Christopher Sims
0019 % Copyright (C) 2010-2011 Dynare Team
0020 %
0021 % This file is part of Dynare.
0022 %
0023 % Dynare is free software: you can redistribute it and/or modify
0024 % it under the terms of the GNU General Public License as published by
0025 % the Free Software Foundation, either version 3 of the License, or
0026 % (at your option) any later version.
0027 %
0028 % Dynare is distributed in the hope that it will be useful,
0029 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0030 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0031 % GNU General Public License for more details.
0032 %
0033 % You should have received a copy of the GNU General Public License
0034 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0035 
0036 
0037 lastwarn('','');
0038 global lq_instruments;
0039 eu=[0;0];C=c;
0040 realsmall=1e-6;
0041 fixdiv=(nargin==6);
0042 n=size(a0,1);
0043 DD=[];E2=[]; E5=0; GAMMA=[];
0044 %
0045 % Find SVD of a0, and create partitions of U, S and V
0046 %
0047 [U0,S0,V0] = svd(a0);
0048 
0049 FL_RANK=rank(S0);
0050 U01=U0(1:n,1:FL_RANK);
0051 U02=U0(1:n,FL_RANK+1:n);
0052 V01=V0(1:n,1:FL_RANK);
0053 V02=V0(1:n,FL_RANK+1:n);
0054 S01=S0(1:FL_RANK,1:FL_RANK);
0055 
0056 C1=U02'*a1*V01;
0057 C2=U02'*a1*V02;
0058 C3=U02'*a2*V01;
0059 C4=U02'*a2*V02;
0060 C5=U02'*PSI;
0061 Sinv=eye(FL_RANK);
0062 for i=1:FL_RANK
0063     Sinv(i,i)=1/S01(i,i);
0064 end
0065 F1=Sinv*U01'*a1*V01;
0066 F2=Sinv*U01'*a1*V02;
0067 F3=Sinv*U01'*a2*V01;
0068 F4=Sinv*U01'*a2*V02;
0069 F5=Sinv*U01'*PSI;
0070 singular=0;
0071 warning('', '');
0072 try
0073     if rcond(C2) < 1e-8
0074         singular=1;
0075     else
0076         warning('off','MATLAB:nearlySingularMatrix');
0077         warning('off','MATLAB:singularMatrix');
0078         UAVinv=inv(C2); % i.e. inv(U02'*a1*V02)
0079         [LastWarningTxt LastWarningID]=lastwarn;
0080         if any(any(isinf(UAVinv)))==1
0081             singular=1;
0082         end
0083     end
0084     if singular == 1 || strcmp('MATLAB:nearlySingularMatrix',LastWarningID) == 1 || ...
0085                  strcmp('MATLAB:illConditionedMatrix',LastWarningID)==1 || ...
0086                  strcmp('MATLAB:singularMatrix',LastWarningID)==1 
0087         [C1,C2,C3,C4, C5, F1, F2, F3, F4, F5, M1, M2, UAVinv, FL_RANK, V01, V02] = PI_gensys_singularC(C1,C2,C3,C4, C5, F1, F2, F3, F4, F5, V01, V02, 0);
0088     end
0089     warning('on','MATLAB:singularMatrix');
0090     warning('on','MATLAB:nearlySingularMatrix');
0091     if (any(any(isinf(UAVinv))) || any(any(isnan(UAVinv)))) 
0092         if(options_.ACES_solver==1)
0093             disp('ERROR! saving PI_gensys_data_dump');
0094             save PI_gensys_data_dump
0095             error('PI_gensys: Inversion of poss. zero matrix UAVinv=inv(U02''*a1*V02)!');
0096         else
0097             warning('PI_gensys: Evading inversion of zero matrix UAVinv=inv(U02''*a1*V02)!');
0098             eu=[0,0];
0099             return;
0100         end
0101     end
0102 catch
0103     errmsg=lasterror;
0104     warning(['error callig PI_gensys_singularC: ' errmsg.message ],'errmsg.identifier');
0105     %error('errcode',['error callig PI_gensys_singularC: ' errmsg.message ]);
0106 end
0107 %
0108 % Define TT1, TT2
0109 %
0110 TT1=V02;
0111 TT2=V01;
0112 
0113 %
0114 %Set up the system matrix for the variables s(t)=V02*Y(t), x(t)=V01*Y(t) and E_t x(t+1)
0115 %                                   and define z(t)'=[s(t)' x(t)]
0116 %
0117 %UAVinv=inv(U02'*a1*V02);
0118 FF=F2; %=Sinv*U01'*a1*V02;
0119 G11=UAVinv*C4;
0120 G12=UAVinv*C3;
0121 G13=-UAVinv*C1;
0122 G31=-FF*G11+F4; % +Sinv*U01'*a2*V02;
0123 G32=-FF*G12+F3; % +Sinv*U01'*a2*V01;
0124 G33=-FF*G13-F1; % -Sinv*U01'*a1*V01;
0125 P1=UAVinv*C5; % *U02'*PSI;
0126 P3=-FF*P1+F5; % + Sinv*U01'*PSI; % This should equal 0
0127 G21=zeros(FL_RANK,(n-FL_RANK));
0128 G22=zeros(FL_RANK,FL_RANK);
0129 G23=eye(FL_RANK);
0130 %H2=zeros(FL_RANK,NX);
0131 num_inst=0;
0132 
0133 % New Definitions
0134 Ze11=zeros(NX,NX); 
0135 Ze12=zeros(NX,(n-FL_RANK)); 
0136 Ze134=zeros(NX,FL_RANK);
0137 Ze31=zeros(FL_RANK,NX);
0138 
0139 % End of New Definitions
0140 
0141 %
0142 % System matrix is called 'G1pi'; Shock matrix is called 'impact'
0143 %
0144 
0145 G1pi=[Ze11 Ze12 Ze134 Ze134; P1 G11 G12 G13; Ze31 G21 G22 G23; P3 G31 G32 G33];
0146 
0147 impact=[eye(NX,NX); zeros(n+FL_RANK,NX)];
0148 
0149 if(options_.ACES_solver==1)
0150     if isfield(lq_instruments,'names')
0151         num_inst=size(lq_instruments.names,1);
0152         if num_inst>0 
0153             i_var=lq_instruments.inst_var_indices;
0154             N1=UAVinv*U02'*lq_instruments.B1;
0155             N3=-FF*N1+Sinv*U01'*lq_instruments.B1;
0156         else
0157             error('WARNING: There are no instrumnets for ACES!');
0158         end
0159         lq_instruments.N1=N1;
0160         lq_instruments.N3=N3;
0161     else
0162         error('WARNING: There are no instrumnets for ACES!');
0163     end
0164     E3=V02*[P1 G11 G12 G13];
0165     E3= E3+ [zeros(size(V01,1),size(E3,2)-size(V01,2)) V01];
0166     E2=-E3;
0167     E5=-V02*N1;
0168     DD=[zeros(NX,size(N1,2));N1; zeros(FL_RANK,size(N1,2));N3];
0169     II=eye(num_inst);
0170     GAMMA=[ E3 -E5 %zeros(size(E3,1),num_inst);
0171             zeros(num_inst,size(E3,2)), II;
0172           ];
0173     eu =[1; 1], nmat=[], gev=[];
0174     return; % do not check B&K compliancy
0175 end
0176 
0177 G0pi=eye(n+FL_RANK+NX);
0178 try
0179     % In Matlab: [aa bb q z v w] = qz(a,b) s.t. qaz = aa, qbz = bb %
0180     % In Octave: [aa bb q z v w] = qz(a,b) s.t. q'az = aa, q'bz=bb %
0181     % and qzcomplex() extension based on lapack zgges produces same
0182     % qz output for Octave as Matlab qz() does for Matlab thus:
0183     if exist('OCTAVE_VERSION')
0184         [a b q z]=qzcomplex(G0pi,G1pi);
0185         q=q';
0186     else
0187         [a b q z]=qz(G0pi,G1pi);
0188     end
0189 catch
0190     try
0191         lerror=lasterror;
0192         disp(['PI_Gensys: ' lerror.message]);
0193         if 0==strcmp('MATLAB:qz:matrixWithNaNInf',lerror.identifier)
0194             disp '** Unexpected Error PI_Gensys:qz: ** :';
0195             button=questdlg('Continue Y/N?','Unexpected Error in qz','No','Yes','Yes'); 
0196             switch button 
0197               case 'No' 
0198                 error ('Terminated')
0199                 %case 'Yes'
0200                 
0201             end
0202         end
0203         G1pi=[];impact=[];nmat=[]; gev=[];
0204         eu=[-2;-2];
0205         return
0206     catch
0207         disp '** Unexpected Error in qz ** :';
0208         disp lerror.message;
0209         button=questdlg('Continue Y/N?','Unexpected Error in qz','No','Yes','Yes'); 
0210         switch button 
0211           case 'No' 
0212             error ('Terminated') 
0213           case 'Yes' 
0214             G1pi=[];impact=[];nmat=[]; gev=[];
0215             eu=[-2;-2];
0216             return
0217         end
0218     end
0219 end
0220 
0221 if ~fixdiv, div=1.01; end
0222 nunstab=0;
0223 zxz=0;
0224 nn=size(a,1);
0225 for i=1:nn
0226     % ------------------div calc------------
0227     if ~fixdiv
0228         if abs(a(i,i)) > 0
0229             divhat=abs(b(i,i))/abs(a(i,i));
0230             % bug detected by Vasco Curdia and Daria Finocchiaro, 2/25/2004  A root of
0231             % exactly 1.01 and no root between 1 and 1.02, led to div being stuck at 1.01
0232             % and the 1.01 root being misclassified as stable.  Changing < to <= below fixes this.
0233             if 1+realsmall<divhat && divhat<=div
0234                 div=.5*(1+divhat);
0235             end
0236         end
0237     end
0238     % ----------------------------------------
0239     nunstab=nunstab+(abs(b(i,i))>div*abs(a(i,i)));
0240     if abs(a(i,i))<realsmall && abs(b(i,i))<realsmall
0241         zxz=1;
0242     end
0243 end
0244 div ;
0245 if ~zxz
0246     [a b q z]=qzdiv(div,a,b,q,z);
0247 end
0248 
0249 gev=[diag(a) diag(b)];
0250 if zxz
0251     disp('Coincident zeros.  Indeterminacy and/or nonexistence.')
0252     eu=[-2;-2];
0253     % correction added 7/29/2003.  Otherwise the failure to set output
0254     % arguments leads to an error message and no output (including eu).
0255     nmat=[]; %;gev=[]
0256     return
0257 end
0258 if (FL_RANK ~= nunstab && options_.ACES_solver~=1)
0259     disp(['Number of unstable variables ' num2str(nunstab)]);
0260     disp( ['does not match number of expectational equations ' num2str(FL_RANK)]); 
0261     nmat=[];% gev=[];
0262     eu=[-2;-2];
0263     return
0264 end
0265 
0266 % New Definitions
0267 z1=z(:,1:n+NX)';
0268 z2=z(:,n+NX+1:n+NX+FL_RANK)';
0269 
0270 % New N Matrix by J Pearlman
0271 z12=z2(:,1:n+NX);
0272 z22=z2(:,n+NX+1:n+NX+FL_RANK);
0273 % End of New Definitions
0274 
0275 % modified by GP:
0276 nmat=real(inv(z22)*z12);
0277 eu=[1;1];