1) relaxed schur_vec_tol

2) warning off divde by zero
3) better figure layout
time-shift
Marco Ratto 2010-02-23 10:47:43 +01:00
parent 263725aa1f
commit bb06ec430a
1 changed files with 82 additions and 40 deletions

View File

@ -21,6 +21,12 @@ function [pdraws, TAU, GAM, LRE, gp, H, JJ] = dynare_identification(options_iden
global M_ options_ oo_ bayestopt_ estim_params_ global M_ options_ oo_ bayestopt_ estim_params_
if exist('OCTAVE_VERSION')
warning('off', 'Octave:divide-by-zero')
else
warning off MATLAB:dividebyzero
end
options_ident = set_default_option(options_ident,'load_ident_files',0); options_ident = set_default_option(options_ident,'load_ident_files',0);
options_ident = set_default_option(options_ident,'useautocorr',1); options_ident = set_default_option(options_ident,'useautocorr',1);
options_ident = set_default_option(options_ident,'ar',3); options_ident = set_default_option(options_ident,'ar',3);
@ -35,7 +41,7 @@ useautocorr = options_ident.useautocorr;
options_.ar=nlags; options_.ar=nlags;
options_.prior_mc = options_ident.prior_mc; options_.prior_mc = options_ident.prior_mc;
options_.options_ident = options_ident; options_.options_ident = options_ident;
options_.Schur_vec_tol = 1.e-8;
options_ = set_default_option(options_,'datafile',[]); options_ = set_default_option(options_,'datafile',[]);
options_.mode_compute = 0; options_.mode_compute = 0;
@ -111,6 +117,11 @@ if iload <=0,
TAU(:,burnin_iteration)=tau; TAU(:,burnin_iteration)=tau;
LRE(:,burnin_iteration)=[oo_.dr.ys(oo_.dr.order_var); vec(g1)]; LRE(:,burnin_iteration)=[oo_.dr.ys(oo_.dr.order_var); vec(g1)];
[gam,stationary_vars] = th_autocovariances(oo0.dr,bayestopt_.mfys,M_,options_); [gam,stationary_vars] = th_autocovariances(oo0.dr,bayestopt_.mfys,M_,options_);
if exist('OCTAVE_VERSION')
warning('off', 'Octave:divide-by-zero')
else
warning off MATLAB:dividebyzero
end
sdy = sqrt(diag(gam{1})); sdy = sqrt(diag(gam{1}));
sy = sdy*sdy'; sy = sdy*sdy';
if useautocorr, if useautocorr,
@ -126,6 +137,7 @@ if iload <=0,
dum = [dum; vec(gam{j+1})]; dum = [dum; vec(gam{j+1})];
end end
GAM(:,burnin_iteration)=[oo_.dr.ys(bayestopt_.mfys); dum]; GAM(:,burnin_iteration)=[oo_.dr.ys(bayestopt_.mfys); dum];
warning warning_old_state;
else else
iteration = iteration + 1; iteration = iteration + 1;
run_index = run_index + 1; run_index = run_index + 1;
@ -206,7 +218,6 @@ if iload <=0,
end end
close(h) close(h)
@ -244,8 +255,11 @@ siLREmean = siLREmean./(max(siLREmean')'*ones(1,estim_params_.np));
tstJmean = derJmean*0; tstJmean = derJmean*0;
tstHmean = derHmean*0; tstHmean = derHmean*0;
tstLREmean = derLREmean*0; tstLREmean = derLREmean*0;
warning_old_state = warning; if exist('OCTAVE_VERSION')
warning off; warning('off', 'Octave:divide-by-zero')
else
warning off MATLAB:dividebyzero
end
for j=1:nparam, for j=1:nparam,
indd = 1:length(siJmean(:,j)); indd = 1:length(siJmean(:,j));
@ -258,7 +272,6 @@ tstLREmean(indd,j-offset) = abs(derLREmean(indd,j-offset))./siLREmean(indd,j-off
end end
end end
warning warning_old_state
if nargout>3 & iload, if nargout>3 & iload,
filnam = dir([IdentifDirectoryName '/' M_.fname '_identif_*.mat']); filnam = dir([IdentifDirectoryName '/' M_.fname '_identif_*.mat']);
@ -373,67 +386,89 @@ disp_identification(pdraws, idemodel, idemoments)
% end % end
% title('Sensitivity in standardized moments'' PCA') % title('Sensitivity in standardized moments'' PCA')
figure, figure('Name','Identification LRE model form'),
subplot(231) subplot(211)
myboxplot(siLREmean) mmm = mean(siLREmean);
[ss, is] = sort(mmm);
myboxplot(siLREmean(:,is))
set(gca,'ylim',[0 1.05]) set(gca,'ylim',[0 1.05])
set(gca,'xticklabel','') set(gca,'xticklabel','')
for ip=1:estim_params_.np, for ip=1:estim_params_.np,
text(ip,-0.02,deblank(M_.param_names(estim_params_.param_vals(ip,1),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') text(ip,-0.02,deblank(M_.param_names(estim_params_.param_vals(is(ip),1),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end end
title('Sensitivity in the LRE model') title('Sensitivity in the LRE model')
subplot(232) subplot(212)
myboxplot(siHmean) mmm = mean(-idelre.Mco');
set(gca,'ylim',[0 1.05]) [ss, is] = sort(mmm);
set(gca,'xticklabel','') myboxplot(idelre.Mco(is,:)')
for ip=1:nparam,
text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Sensitivity in the model')
subplot(233)
myboxplot(siJmean)
set(gca,'ylim',[0 1.05])
set(gca,'xticklabel','')
for ip=1:nparam,
text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Sensitivity in the moments')
subplot(234)
myboxplot(idelre.Mco')
set(gca,'ylim',[0 1]) set(gca,'ylim',[0 1])
set(gca,'xticklabel','') set(gca,'xticklabel','')
for ip=1:estim_params_.np, for ip=1:estim_params_.np,
text(ip,-0.02,deblank(M_.param_names(estim_params_.param_vals(ip,1),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') text(ip,-0.02,deblank(M_.param_names(estim_params_.param_vals(is(ip),1),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end end
title('Multicollinearity in the LRE model') title('Multicollinearity in the LRE model')
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_LRE'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_LRE']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_LRE']);
if options_.nograph, close(gcf); end
subplot(235) figure('Name','Identification in the model'),
myboxplot(idemodel.Mco') subplot(211)
mmm = mean(siHmean);
[ss, is] = sort(mmm);
myboxplot(siHmean(:,is))
set(gca,'ylim',[0 1.05])
set(gca,'xticklabel','')
for ip=1:nparam,
text(ip,-0.02,bayestopt_.name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Sensitivity in the model')
subplot(212)
mmm = mean(-idemodel.Mco');
[ss, is] = sort(mmm);
myboxplot(idemodel.Mco(is,:)')
set(gca,'ylim',[0 1]) set(gca,'ylim',[0 1])
set(gca,'xticklabel','') set(gca,'xticklabel','')
for ip=1:nparam, for ip=1:nparam,
text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none') text(ip,-0.02,bayestopt_.name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end end
title('Multicollinearity in the model') title('Multicollinearity in the model')
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_model'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_model']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_model']);
if options_.nograph, close(gcf); end
subplot(236) figure('Name','Identification in the moments'),
myboxplot(idemoments.Mco') subplot(211)
mmm = mean(siJmean);
[ss, is] = sort(mmm);
myboxplot(siJmean(:,is))
set(gca,'ylim',[0 1.05])
set(gca,'xticklabel','')
for ip=1:nparam,
text(ip,-0.02,bayestopt_.name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Sensitivity in the moments')
subplot(212)
mmm = mean(-idemoments.Mco');
[ss, is] = sort(mmm);
myboxplot(idemoments.Mco(is,:)')
set(gca,'ylim',[0 1]) set(gca,'ylim',[0 1])
set(gca,'xticklabel','') set(gca,'xticklabel','')
for ip=1:nparam, for ip=1:nparam,
text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none') text(ip,-0.02,bayestopt_.name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end end
title('Multicollinearity in the moments') title('Multicollinearity in the moments')
saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident']) saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_moments'])
eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident']); eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_moments']);
eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident']); eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_moments']);
if options_.nograph, close(gcf); end if options_.nograph, close(gcf); end
figure, figure('Name','Condition Number'),
subplot(221) subplot(221)
hist(log10(idemodel.cond)) hist(log10(idemodel.cond))
title('log10 of Condition number in the model') title('log10 of Condition number in the model')
@ -527,3 +562,10 @@ for j=1:nparam,
if options_.nograph, close(gcf); end if options_.nograph, close(gcf); end
end end
end end
if exist('OCTAVE_VERSION')
warning('on', 'Octave:divide-by-zero')
else
warning on MATLAB:dividebyzero
end