diff --git a/matlab/partial_information/PI_gensys.m b/matlab/partial_information/PI_gensys.m
index 6cb64f702..fd4d8255b 100644
--- a/matlab/partial_information/PI_gensys.m
+++ b/matlab/partial_information/PI_gensys.m
@@ -34,7 +34,7 @@ function [G1pi,C,impact,nmat,TT1,TT2,gev,eu, DD, E2, E5, GAMMA, FL_RANK ]=PI_gen
% along with Dynare. If not, see .
-lastwarn('');
+lastwarn('','');
global lq_instruments;
eu=[0;0];C=c;
realsmall=1e-6;
@@ -67,17 +67,23 @@ 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
+singular=0;
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);
+ if rcond(C2) < 1e-8
+ singular=1;
+ else
+ warning('off','MATLAB:nearlySingularMatrix');
+ warning('off','MATLAB:singularMatrix');
+ UAVinv=inv(C2); % i.e. inv(U02'*a1*V02)
+ [LastWarningTxt LastWarningID]=lastwarn;
+ if any(any(isinf(UAVinv)))==1
+ singular=1;
+ end
+ end
+ if singular == 1 || strcmp('MATLAB:nearlySingularMatrix',LastWarningID) == 1 || ...
+ strcmp('MATLAB:illConditionedMatrix',LastWarningID)==1 || ...
+ strcmp('MATLAB:singularMatrix',LastWarningID)==1
[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;
@@ -95,11 +101,10 @@ try
return;
end
end
-%
catch
[errmsg, errcode]=lasterror;
warning(['error callig PI_gensys_singularC: ' errmsg ],'errcode');
- warning('','');
+ error('errcode',['error callig PI_gensys_singularC: ' errmsg ]);
end
%
% Define TT1, TT2
diff --git a/matlab/partial_information/PI_gensys_singularC.m b/matlab/partial_information/PI_gensys_singularC.m
index 1de00c30c..d95232540 100644
--- a/matlab/partial_information/PI_gensys_singularC.m
+++ b/matlab/partial_information/PI_gensys_singularC.m
@@ -24,7 +24,7 @@ function [C1,C2,C3,C4, C5, F1, F2, F3, F4, F5, M1, M2, UAVinv, FL_RANK]=PI_gensy
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see .
-level=level+1;
+level=level+1
if level>100
error( ' PI_gensys_singularC recurssion exceeeded its maximum of 100 iterations! ');
end
@@ -61,18 +61,27 @@ F3 = M*F3*M1;
F2 =[M*F2 M*F1*M2];
F1 = M*F1*M1;
warning('', '');
+singular=0;
try
- UAVinv=inv(C2);
- [LastWarningTxt LastWarningID]=lastwarn;
- if strcmp('MATLAB:nearlySingularMatrix',LastWarningID) | ...
- strcmp('MATLAB:illConditionedMatrix',LastWarningID) | ...
- strcmp('MATLAB:singularMatrix',LastWarningID) | isinf(UAVinv)
+ if rcond(C2) < 1e-8
+ singular=1;
+ else
+ UAVinv=inv(C2);
+ [LastWarningTxt LastWarningID]=lastwarn;
+ if any(any(isinf(UAVinv)))==1
+ singular=1;
+ end
+ end
+ % line test is for Octave strncmp('warning: inverse: matrix singular',LastWarningTxt, 33)==1 || ...
+ if singular==1 || strcmp('MATLAB:nearlySingularMatrix',LastWarningID)==1 || ...
+ strcmp('MATLAB:illConditionedMatrix',LastWarningID)==1 || ...
+ strcmp('MATLAB:singularMatrix',LastWarningID)==1
[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;
+ [errmsg, errcode]=lasterr;
warning(['error callig PI_gensys_singularC: ' errmsg ],'errcode');
- warning('','');
+ error('errcode',['error callig PI_gensys_singularC: ' errmsg ]);
end
return;
diff --git a/matlab/partial_information/dr1_PI.m b/matlab/partial_information/dr1_PI.m
index aa931b397..fbab7842b 100644
--- a/matlab/partial_information/dr1_PI.m
+++ b/matlab/partial_information/dr1_PI.m
@@ -109,7 +109,7 @@ function [dr,info,M_,options_,oo_] = dr1_PI(dr,task,M_,options_,oo_)
end
- if options_.ACES_solver == 1
+ if options_.ACES_solver == 1
sim_ruleids=[];
tct_ruleids=[];
if size(M_.equations_tags,1)>0 % there are tagged equations, check if they are aceslq rules
@@ -124,37 +124,37 @@ function [dr,info,M_,options_,oo_] = dr1_PI(dr,task,M_,options_,oo_)
end
lq_instruments.sim_ruleids=sim_ruleids;
lq_instruments.tct_ruleids=tct_ruleids;
- %if isfield(lq_instruments,'xsopt_SS') %% changed by BY
+ %if isfield(lq_instruments,'xsopt_SS') %% changed by BY
[junk, lq_instruments.xsopt_SS,lq_instruments.lmopt_SS,s2,check] = opt_steady_get;%% changed by BY
[qc, DYN_Q] = QPsolve(lq_instruments, s2, check); %% added by BY
z = repmat(lq_instruments.xsopt_SS,1,klen);
- else
+ else
z = repmat(dr.ys,1,klen);
end
z = z(iyr0) ;
[junk,jacobia_] = feval([M_.fname '_dynamic'],z,[oo_.exo_simul ...
- oo_.exo_det_simul], M_.params, it_);
-
- if options_.ACES_solver==1 & (length(sim_ruleids)>0 || length(tct_ruleids)>0 )
- if length(sim_ruleids)>0
- sim_rule=jacobia_(sim_ruleids,:);
- % uses the subdirectory - BY
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_sim_rule.txt'], 'sim_rule', '-ascii', '-double', '-tabs');
- end
- if length(tct_ruleids)>0
- tct_rule=jacobia_(tct_ruleids,:);
- % uses the subdirectory - BY
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_tct_rule.txt'], 'tct_rule', '-ascii', '-double', '-tabs');
- end
- aces_ruleids=union(tct_ruleids,sim_ruleids);
- j_size=size(jacobia_,1);
- j_rows=1:j_size;
- j_rows = setxor(j_rows,aces_ruleids);
- jacobia_=jacobia_(j_rows ,:);
- end
-
+ oo_.exo_det_simul], M_.params, it_);
+
+ if options_.ACES_solver==1 && (length(sim_ruleids)>0 || length(tct_ruleids)>0 )
+ if length(sim_ruleids)>0
+ sim_rule=jacobia_(sim_ruleids,:);
+ % uses the subdirectory - BY
+ save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_sim_rule.txt'], 'sim_rule', '-ascii', '-double', '-tabs');
+ end
+ if length(tct_ruleids)>0
+ tct_rule=jacobia_(tct_ruleids,:);
+ % uses the subdirectory - BY
+ save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_tct_rule.txt'], 'tct_rule', '-ascii', '-double', '-tabs');
+ end
+ aces_ruleids=union(tct_ruleids,sim_ruleids);
+ j_size=size(jacobia_,1);
+ j_rows=1:j_size;
+ j_rows = setxor(j_rows,aces_ruleids);
+ jacobia_=jacobia_(j_rows ,:);
+ end
+
end
-
+
if options_.debug
save([M_.fname '_debug.mat'],'jacobia_')
end