Partial Information: Update with improved, generalised recursive PCL solution method

matlab/partial_information/PI_gensys.m

-function [G1pi,C,impact,nmat,TT1,TT2,gev,eu, DD, E3, E5]=PI_gensys(a0,a1,a2,a3,c,PSI,NX,NETA,FL_RANK,M_,options_)
+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_)
 % 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.
@@ -34,12 +34,13 @@ function [G1pi,C,impact,nmat,TT1,TT2,gev,eu, DD, E3, E5]=PI_gensys(a0,a1,a2,a3,c
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 
+lastwarn('');
 global lq_instruments;
 eu=[0;0];C=c;
 realsmall=1e-6;
 fixdiv=(nargin==6);
 n=size(a0,1);
-DD=[];E3=[]; E5=0;%[];
+DD=[];E2=[]; E5=0; GAMMA=[];
 %
 % Find SVD of a0, and create partitions of U, S and V
 %
@@ -51,84 +52,84 @@ U02=U0(1:n,FL_RANK+1:n);
 V01=V0(1:n,1:FL_RANK);
 V02=V0(1:n,FL_RANK+1:n);
 S01=S0(1:FL_RANK,1:FL_RANK);
+
+C1=U02'*a1*V01;
+C2=U02'*a1*V02;
+C3=U02'*a2*V01;
+C4=U02'*a2*V02;
+C5=U02'*PSI;
+Sinv=eye(FL_RANK);
+for i=1:FL_RANK
+    Sinv(i,i)=1/S01(i,i);
+end
+F1=Sinv*U01'*a1*V01;
+F2=Sinv*U01'*a1*V02;
+F3=Sinv*U01'*a2*V01;
+F4=Sinv*U01'*a2*V02;
+F5=Sinv*U01'*PSI;
+%if isempty(V02) & isempty(U02) % no only backward looking variables model
+warning('', '');
+try
+  warning('off','MATLAB:nearlySingularMatrix');
+  warning('off','MATLAB:singularMatrix');
+  UAVinv=inv(C2); % i.e. inv(U02'*a1*V02)
+  [LastWarningTxt LastWarningID]=lastwarn;
+  if (strcmp('MATLAB:nearlySingularMatrix',LastWarningID) | ...
+      strcmp('MATLAB:illConditionedMatrix',LastWarningID) | ...
+      strcmp('MATLAB:singularMatrix',LastWarningID) | isinf(UAVinv))
+    %display(LastWarningTxt);
+    [C1,C2,C3,C4, C5, F1, F2, F3, F4, F5, M1, M2, UAVinv, FL_RANK] = PI_gensys_singularC(C1,C2,C3,C4, C5, F1, F2, F3, F4, F5, 0);
+    V02=[V02 V01*M2];
+    V01=V01*M1;
+  end
+  warning('on','MATLAB:singularMatrix');
+  warning('on','MATLAB:nearlySingularMatrix');
+  if (isinf(UAVinv) | isnan(UAVinv)) 
+    if(options_.useACES==1)
+      disp('ERROR! saving PI_gensys_data_dump');
+      save PI_gensys_data_dump
+      error('PI_gensys: Inversion of poss. zero matrix UAVinv=inv(U02''*a1*V02)!');
+    else
+      warning('PI_gensys: Evading inversion of zero matrix UAVinv=inv(U02''*a1*V02)!');
+      eu=[0,0];
+      return;
+    end
+  end
+%
+catch
+  [errmsg, errcode]=lasterror;
+  warning(['error callig PI_gensys_singularC: ' errmsg ],'errcode');
+  warning('','');
+end
 %
 % Define TT1, TT2
 %
 TT1=V02;
 TT2=V01;
-%
-%Invert S01
-%
-Sinv=eye(FL_RANK);
-for i=1:FL_RANK
-    Sinv(i,i)=1/S01(i,i);
-end
+
 %
 %Set up the system matrix for the variables s(t)=V02*Y(t), x(t)=V01*Y(t) and E_t x(t+1)
 %                                   and define z(t)'=[s(t)' x(t)]
 %
-FF=Sinv*U01'*a1*V02;
-UAVinv=inv(U02'*a1*V02);
-G11=UAVinv*U02'*a2*V02;
-G12=UAVinv*U02'*a2*V01;
-G13=-UAVinv*U02'*a1*V01;
-G31=-FF*G11+Sinv*U01'*a2*V02;
-G32=-FF*G12+Sinv*U01'*a2*V01;
-G33=-FF*G13-Sinv*U01'*a1*V01;
-H1=UAVinv*U02'*PSI;
-H3=-FF*H1+Sinv*U01'*PSI; % This should equal 0
-G21=zeros(FL_RANK,n-FL_RANK);
+%UAVinv=inv(U02'*a1*V02);
+FF=F2; %=Sinv*U01'*a1*V02;
+G11=UAVinv*C4;
+G12=UAVinv*C3;
+G13=-UAVinv*C1;
+G31=-FF*G11+F4; % +Sinv*U01'*a2*V02;
+G32=-FF*G12+F3; % +Sinv*U01'*a2*V01;
+G33=-FF*G13-F1; % -Sinv*U01'*a1*V01;
+P1=UAVinv*C5; % *U02'*PSI;
+P3=-FF*P1+F5; % + Sinv*U01'*PSI; % This should equal 0
+G21=zeros(FL_RANK,(n-FL_RANK));
 G22=zeros(FL_RANK,FL_RANK);
 G23=eye(FL_RANK);
 %H2=zeros(FL_RANK,NX);
 num_inst=0;
-P1=H1;
-P3=H3;
-if(options_.ACES_solver==1 & isfield(lq_instruments,'names'))
-    num_inst=size(lq_instruments.names,1);
-    if ~isfield(lq_instruments,'inst_varexo_indices') & num_inst>0
-      for i=1:num_inst
-        i_tmp = strmatch(deblank(lq_instruments.names(i,:)),M_.exo_names,'exact');
-        if isempty(i_tmp)
-          error (['One of the specified instrument variables does not exist']) ;
-        else
-          i_var(i) = i_tmp;
-        end
-      end
-      lq_instruments.inst_varexo_indices=i_var;
-    elseif size(lq_instruments.inst_varexo_indices)>0
-        i_var=lq_instruments.inst_varexo_indices;
-        if ~num_inst
-            num_inst=size(lq_instruments.inst_varexo_indices);
-        end
-    else
-        i_var=[];
-        num_inst=0;
-    end
-    if size(i_var,2)>0 & size(i_var,2)==num_inst
-        N1=H1(:,i_var);
-        N3=H3(:,i_var);
-        x_var=zeros(NX,1);
-        for i=1:NX
-            if isempty(find(i_var==i))
-                x_var(i)=i;
-            end
-        end
-        x_var=nonzeros(x_var);
-        P1=H1(:,x_var);
-        P3=H3(:,x_var);
-        NX=NX-num_inst;
-    else
-        error('WARNING: There are no instrumnets for ACES!');
-    end
-    lq_instruments.N1=N1;
-    lq_instruments.N3=N3;
-elseif(options_.ACES_solver==1)
-    error('WARNING: There are no instrumnets for ACES!');
-end
+
 % New Definitions
 Ze11=zeros(NX,NX); 
-Ze12=zeros(NX,n-FL_RANK); 
+Ze12=zeros(NX,(n-FL_RANK)); 
 Ze134=zeros(NX,FL_RANK);
 Ze31=zeros(FL_RANK,NX);
 
@@ -143,12 +144,31 @@ G1pi=[Ze11 Ze12 Ze134 Ze134; P1 G11 G12 G13; Ze31 G21 G22 G23; P3 G31 G32 G33];
 impact=[eye(NX,NX); zeros(n+FL_RANK,NX)];
 
 if(options_.ACES_solver==1)
-    E3=V02*[P1 G11 G12 G13];
-    E3=E3+ [zeros(size(V01,1),size(E3,2)-size(V01,2)) V01];
-    E5=-V02*N1;
-    DD=[zeros(NX,size(N1,2));N1; zeros(FL_RANK,size(N1,2));N3];
-    eu =[1; 1], nmat=[], gev=[];
-    return; % do not check B&K compliancy
+  if isfield(lq_instruments,'names')
+    num_inst=size(lq_instruments.names,1);
+    if num_inst>0 
+      i_var=lq_instruments.inst_var_indices;
+      N1=UAVinv*U02'*lq_instruments.B1;
+      N3=-FF*N1+Sinv*U01'*lq_instruments.B1;
+    else
+        error('WARNING: There are no instrumnets for ACES!');
+    end
+    lq_instruments.N1=N1;
+    lq_instruments.N3=N3;
+  else
+    error('WARNING: There are no instrumnets for ACES!');
+  end
+  E3=V02*[P1 G11 G12 G13];
+  E3= E3+ [zeros(size(V01,1),size(E3,2)-size(V01,2)) V01];
+  E2=-E3;
+  E5=-V02*N1;
+  DD=[zeros(NX,size(N1,2));N1; zeros(FL_RANK,size(N1,2));N3];
+  II=eye(num_inst);
+  GAMMA=[ E3 -E5 %zeros(size(E3,1),num_inst);
+    zeros(num_inst,size(E3,2)), II;
+    ];
+  eu =[1; 1], nmat=[], gev=[];
+  return; % do not check B&K compliancy
 end
 
 G0pi=eye(n+FL_RANK+NX);

matlab/partial_information/PI_gensys_singularC.m

+function [C1,C2,C3,C4, C5, F1, F2, F3, F4, F5, M1, M2, UAVinv, FL_RANK]=PI_gensys_singularC(C1in, C2in, C3in, C4in, C5in, F1, F2, F3, F4, F5, level)
+% [C1,C2,C3,C4, C5, F1, F2, F3, F4, F5, M1, M2, UAVinv,FL_RANK]...
+%         =PI_gensys_singularC(C1in, C2in, C3in, C4in, C5in, F1, F2, F3, F4, F5, level)
+% 
+% Recursive extension for PI_gensys function PCL general DSGE solver
+% devised by Prof. Joseph Pearlman 
+% developed by George Perendia
+% December 2010
+
+% Copyright (C) 1996-2011 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+level=level+1;
+if level>100
+  error( ' PI_gensys_singularC recurssion exceeeded its maximum of 100 iterations! ');
+end
+warning('', '');
+M1=[];M2=[]; UAVinv=[];
+%
+% Find SVD of a0, and create partitions of U, S and V
+%
+
+[J0,K0,L0] = svd(C2in);
+n=size(C2in,1);
+K_RANK=rank(K0);
+J2=J0(1:n,K_RANK+1:n);
+
+J2C1=J2'*C1in;
+M = null(J2C1)';
+MJCinv= inv([M;J2C1]);
+[sm1, sm2]=size (M');
+M1=MJCinv(1:sm1,1:sm2);
+M2=MJCinv(1:sm1,1+sm2:end);
+FL_RANK=rank(M);
+
+%Define new Cs
+C5=[ C5in; J2C1*F5];
+C4=[C4in C3in*M2; J2C1*F4 J2C1*F3*M2];
+C3=[ C3in*M1; J2C1*F3*M1];
+C2=[C2in C1in*M2; J2'*(C4in+C1in*F2) J2'*(C3in+C1in*F1)*M2];
+C1=[ C1in*M1; J2'*(C3in+C1in*F1)*M1];
+%define new after Cs Fs
+% keep this reverse order!!
+F5 =M*F5 ;
+F4 =[M*F4  M*F3*M2];
+F3 = M*F3*M1;
+F2 =[M*F2  M*F1*M2];
+F1 = M*F1*M1;
+warning('', '');
+try
+  UAVinv=inv(C2);
+  [LastWarningTxt LastWarningID]=lastwarn;
+  if strcmp('MATLAB:nearlySingularMatrix',LastWarningID) | ...
+      strcmp('MATLAB:illConditionedMatrix',LastWarningID) | ...
+      strcmp('MATLAB:singularMatrix',LastWarningID) | isinf(UAVinv)
+    [C1,C2,C3,C4, C5, F1, F2, F3, F4, F5, M1, M2, UAVinv, FL_RANK] = PI_gensys_singularC(C1,C2,C3,C4, C5, F1, F2, F3, F4, F5, level);
+  end
+catch
+  [errmsg, errcode]=lasterror;
+  warning(['error callig PI_gensys_singularC: ' errmsg ],'errcode');
+  warning('','');
+end
+
+return;
+