stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Partial information changes that make adjustment for Octave and use rcond() < 1e-8 ..." to determine if a matrix is invertible plus some minor bug and formatting changes in dr1_PI.m

time-shift
George Perendia 2011-01-17 20:38:49 +00:00
parent 6b1f84ac27
commit 595675d333
3 changed files with 58 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
matlab/partial_information/PI_gensys.m
View File

@ -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 <http://www.gnu.org/licenses/>.
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

25
matlab/partial_information/PI_gensys_singularC.m
View File

@ -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 <http://www.gnu.org/licenses/>.
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;

48
matlab/partial_information/dr1_PI.m
View File

@ -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