From 57eb9e12e5c82d8dba5b09254256f5bbf1e7c82d Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Tue, 10 Dec 2013 12:42:45 +0100 Subject: [PATCH 1/6] 1) Allow for endogenous prior restrictions 2) trap cases where only very few parameter combinations provide unique saddle path solution --- matlab/gsa/stab_map_.m | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/matlab/gsa/stab_map_.m b/matlab/gsa/stab_map_.m index f2d47c78e..44b006a1f 100644 --- a/matlab/gsa/stab_map_.m +++ b/matlab/gsa/stab_map_.m @@ -257,6 +257,9 @@ if fload==0, %try stoch_simul([]); try [Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict'); + if info(1)==0, + info=endogenous_prior(Tt,Rr,M_,options_,oo_); + end infox(j,1)=info(1); if infox(j,1)==0 && ~exist('T'), dr_=oo_.dr; @@ -533,6 +536,9 @@ if length(iunstable)>0 && length(iunstable)0 && length(iunstable)10, + stab_map_2(lpmat(istable,:),alpha2, pvalue_corr, asname, OutputDirectoryName,xparam1,astitle); + end if length(iunstable)>10, stab_map_2(lpmat(iunstable,:),alpha2, pvalue_corr, auname, OutputDirectoryName,xparam1,autitle); end From 5516c6420511f0b33791348074bffd3d6d352202 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Tue, 10 Dec 2013 14:10:43 +0100 Subject: [PATCH 2/6] Fixed bug introduced in commit 57eb9e12e5c82d8dba5b09254256f5bbf1e7c82d --- matlab/gsa/stab_map_.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/gsa/stab_map_.m b/matlab/gsa/stab_map_.m index 44b006a1f..a847c68ef 100644 --- a/matlab/gsa/stab_map_.m +++ b/matlab/gsa/stab_map_.m @@ -258,7 +258,7 @@ if fload==0, try [Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict'); if info(1)==0, - info=endogenous_prior(Tt,Rr,M_,options_,oo_); + info=endogenous_prior_restrictions(Tt,Rr,M_,options_,oo_); end infox(j,1)=info(1); if infox(j,1)==0 && ~exist('T'), From f8ec220711c0bf173796d2188807386ba109a967 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Tue, 10 Dec 2013 14:32:09 +0100 Subject: [PATCH 3/6] 1) Fixed bug introduced in commit 3cb16af5e323ada152cbac0a2a55cd4607bf6e7b 2) Harmonize bivariate plots in the proper way 3) Harmonize bivariate plots for rmse analysis --- matlab/gsa/filt_mc_.m | 22 +++++++++++----------- matlab/gsa/stab_map_.m | 20 ++++++++++---------- matlab/gsa/stab_map_2.m | 7 +------ 3 files changed, 22 insertions(+), 27 deletions(-) diff --git a/matlab/gsa/filt_mc_.m b/matlab/gsa/filt_mc_.m index b4bc83af5..ed6729362 100644 --- a/matlab/gsa/filt_mc_.m +++ b/matlab/gsa/filt_mc_.m @@ -524,17 +524,17 @@ else % plot trade-offs if ~options_.nograph a00=jet(size(vvarvecm,1)); - for ix=1:ceil(length(nsnam)/5), - if options_.opt_gsa.ppost - temp_name='RMSE Posterior Tradeoffs: Log Posterior'; + if options_.opt_gsa.ppost + temp_name='RMSE Posterior Tradeoffs:'; + else + if options_.opt_gsa.pprior + temp_name='RMSE Prior Tradeoffs:'; else - if options_.opt_gsa.pprior - temp_name='RMSE Prior Tradeoffs: Log Posterior'; - else - temp_name='RMSE MC Tradeoffs: Log Posterior'; - end - end - hh = dyn_figure(options_,'name',[temp_name,' ',int2str(ix)]); + temp_name='RMSE MC Tradeoffs;'; + end + end + for ix=1:ceil(length(nsnam)/5), + hh = dyn_figure(options_,'name',[temp_name,' observed variables ',int2str(ix)]); for j=1+5*(ix-1):min(size(snam2,1),5*ix), subplot(2,3,j-5*(ix-1)) %h0=cumplot(x(:,nsnam(j)+nshock)); @@ -617,7 +617,7 @@ else fnam = ['rmse_mc_',deblank(vvarvecm(i,:))]; end end - stab_map_2(x(ixx(1:nfilt0(i),i),:),alpha2,pvalue,fnam, OutDir); + stab_map_2(x(ixx(1:nfilt0(i),i),:),alpha2,pvalue,fnam, OutDir,[],[temp_name ' observed variable ' deblank(vvarvecm(i,:))]); % [pc,latent,explained] = pcacov(c0); % %figure, bar([explained cumsum(explained)]) diff --git a/matlab/gsa/stab_map_.m b/matlab/gsa/stab_map_.m index a847c68ef..c184c20aa 100644 --- a/matlab/gsa/stab_map_.m +++ b/matlab/gsa/stab_map_.m @@ -471,11 +471,11 @@ if pprior aunstablename=[aname, '_unst']; aunstabletitle='Prior StabMap: Parameter driving explosiveness of solution'; awronguniname=[aname, '_wrong']; awrongunititle='Prior StabMap: Parameter driving inability to find solution'; % bivariate - auname='prior_unacceptable'; autitle='Prior non-existence of unique stable solution (Unacceptable)'; - aunstname='prior_unstable'; aunsttitle='Prior explosiveness of solution'; - aindname='prior_indeterm'; aindtitle='Prior Indeterminacy'; - awrongname='prior_wrong'; awrongtitle='Prior inability to find solution'; - asname='prior_stable'; astitle='Prior unique Stable Saddle-Path'; + auname='prior_unacceptable'; autitle='Prior StabMap: non-existence of unique stable solution (Unacceptable)'; + aunstname='prior_unstable'; aunsttitle='Prior StabMap: explosiveness of solution'; + aindname='prior_indeterm'; aindtitle='Prior StabMap: Indeterminacy'; + awrongname='prior_wrong'; awrongtitle='Prior StabMap: inability to find solution'; + asname='prior_stable'; astitle='Prior StabMap: unique Stable Saddle-Path'; else % univariate aname='mc_stab'; atitle='MC (around posterior mode) StabMap: Parameter driving non-existence of unique stable solution (Unacceptable)'; @@ -483,11 +483,11 @@ else aunstablename=[aname, '_unst']; aunstabletitle='MC (around posterior mode) StabMap: Parameter driving explosiveness of solution'; awronguniname=[aname, '_wrong']; awrongunititle='MC (around posterior mode) StabMap: Parameter driving inability to find solution'; % bivariate - auname='mc_unacceptable'; autitle='MC (around posterior mode) non-existence of unique stable solution (Unacceptable)'; - aunstname='mc_unstable'; aunsttitle='MC (around posterior mode) explosiveness of solution'; - aindname='mc_indeterm'; aindtitle='MC (around posterior mode) Indeterminacy'; - awrongname='mc_wrong'; awrongtitle='MC (around posterior mode) inability to find solution'; - asname='mc_stable'; astitle='MC (around posterior mode) Unique Stable Saddle-Path'; + auname='mc_unacceptable'; autitle='MC (around posterior mode) StabMap: non-existence of unique stable solution (Unacceptable)'; + aunstname='mc_unstable'; aunsttitle='MC (around posterior mode) StabMap: explosiveness of solution'; + aindname='mc_indeterm'; aindtitle='MC (around posterior mode) StabMap: Indeterminacy'; + awrongname='mc_wrong'; awrongtitle='MC (around posterior mode) StabMap: inability to find solution'; + asname='mc_stable'; astitle='MC (around posterior mode) StabMap: Unique Stable Saddle-Path'; end delete([OutputDirectoryName,filesep,fname_,'_',aname,'_*.*']); %delete([OutputDirectoryName,filesep,fname_,'_',aname,'_SA_*.*']); diff --git a/matlab/gsa/stab_map_2.m b/matlab/gsa/stab_map_2.m index e997af919..d4e94f02e 100644 --- a/matlab/gsa/stab_map_2.m +++ b/matlab/gsa/stab_map_2.m @@ -86,14 +86,9 @@ for j=1:npar, fprintf(1,'%20s: corrcoef = %7.3f\n',tmp_name,c0(i2(jx),j)); if ~options_.nograph, - if strcmp(fnam(1:2),'mc') - type='MC (around posterior mode) StabMap: '; - elseif strcmp(fnam(1:5),'prior') - type='Prior StabMap: '; - end if mod(j2,12)==1, ifig=ifig+1; - hh=dyn_figure(options_,'name',[type,'Correlations in the ',figtitle,' sample ', num2str(ifig)]); + hh=dyn_figure(options_,'name',[figtitle,' sample bivariate projection ', num2str(ifig)]); end subplot(3,4,j2-(ifig-1)*12) % bar(c0(i2,j)), From db9b61df7826fa8666f5a2831e379a535b929105 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Tue, 10 Dec 2013 17:10:04 +0100 Subject: [PATCH 4/6] Bug fixes, reduced the number of graphs and improved info on graphical output for screening identification and MC theoretical moments in GSA. --- matlab/gsa/map_ident_.m | 68 ++++++++++++++++++++--------------------- matlab/gsa/mc_moments.m | 12 +++++--- matlab/gsa/th_moments.m | 7 +++-- 3 files changed, 46 insertions(+), 41 deletions(-) diff --git a/matlab/gsa/map_ident_.m b/matlab/gsa/map_ident_.m index 591c547c2..867dfcfd6 100644 --- a/matlab/gsa/map_ident_.m +++ b/matlab/gsa/map_ident_.m @@ -209,7 +209,7 @@ if opt_gsa.morris==1, load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAvdec','vdec','ir_vdec','ic_vdec') end - hh = dyn_figure(options_); + hh = dyn_figure(options_,'name','Screening identification: variance decomposition'); % boxplot(SAvdec,'whis',10,'symbol','r.') myboxplot(SAvdec,[],'.',[],10) set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT]) @@ -221,7 +221,7 @@ if opt_gsa.morris==1, text(ip,-2,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none') end xlabel(' ') - title('All variance decomposition') + title('Elementary effects variance decomposition') dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morris_vdec'],options_); else save([OutputDirectoryName,'/',fname_,'_morris_IDE'],'vdec') @@ -314,7 +314,7 @@ if opt_gsa.morris==1, load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'ac','ir_ac','ic_ac') end - hh=dyn_figure(options_); + hh=dyn_figure(options_,'name','Screening identification: theoretical moments'); % boxplot(SAcc,'whis',10,'symbol','r.') myboxplot(SAcc,[],'.',[],10) set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT]) @@ -326,7 +326,7 @@ if opt_gsa.morris==1, text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none') end xlabel(' ') - title('EET All moments') + title('Elementary effects in the moments') dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morris_moments'],options_); % close(gcf), @@ -709,7 +709,7 @@ if opt_gsa.morris==1, else load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAnorm','SAmunorm','SAsignorm') end - hh=dyn_figure(options_); %bar(SAnorm(:,irel)) + hh=dyn_figure(options_,'name','Screening identification: model'); %bar(SAnorm(:,irel)) % boxplot(SAnorm','whis',10,'symbol','r.') myboxplot(SAnorm',[],'.',[],10) set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT]) @@ -725,35 +725,35 @@ if opt_gsa.morris==1, title('Elementary effects in the model') dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morris_par'],options_); - hh=dyn_figure(options_); %bar(SAmunorm(:,irel)) -% boxplot(SAmunorm','whis',10,'symbol','r.') - myboxplot(SAmunorm',[],'.',[],10) - set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT]) - set(gca,'xlim',[0.5 npT+0.5]) - set(gca,'ylim',[-1 1]) - set(gca,'position',[0.13 0.2 0.775 0.7]) - xlabel(' ') - for ip=1:npT, - text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none') - end - xlabel(' ') - title('\mu in the model') - dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morrismu_par'],options_); - - hh=dyn_figure(options_); %bar(SAsignorm(:,irel)) -% boxplot(SAsignorm','whis',10,'symbol','r.') - myboxplot(SAsignorm',[],'.',[],10) - set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT]) - set(gca,'xlim',[0.5 npT+0.5]) - set(gca,'ylim',[0 1]) - set(gca,'position',[0.13 0.2 0.775 0.7]) - xlabel(' ') - for ip=1:npT, - text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none') - end - xlabel(' ') - title('\sigma in the model') - dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morrissig_par'],options_); +% hh=dyn_figure(options_); %bar(SAmunorm(:,irel)) +% % boxplot(SAmunorm','whis',10,'symbol','r.') +% myboxplot(SAmunorm',[],'.',[],10) +% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT]) +% set(gca,'xlim',[0.5 npT+0.5]) +% set(gca,'ylim',[-1 1]) +% set(gca,'position',[0.13 0.2 0.775 0.7]) +% xlabel(' ') +% for ip=1:npT, +% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none') +% end +% xlabel(' ') +% title('\mu in the model') +% dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morrismu_par'],options_); +% +% hh=dyn_figure(options_); %bar(SAsignorm(:,irel)) +% % boxplot(SAsignorm','whis',10,'symbol','r.') +% myboxplot(SAsignorm',[],'.',[],10) +% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT]) +% set(gca,'xlim',[0.5 npT+0.5]) +% set(gca,'ylim',[0 1]) +% set(gca,'position',[0.13 0.2 0.775 0.7]) +% xlabel(' ') +% for ip=1:npT, +% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none') +% end +% xlabel(' ') +% title('\sigma in the model') +% dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morrissig_par'],options_); % figure, bar(SAnorm(:,irel)') % set(gca,'xtick',[1:j0]) diff --git a/matlab/gsa/mc_moments.m b/matlab/gsa/mc_moments.m index 31fd3bce6..4b0624145 100644 --- a/matlab/gsa/mc_moments.m +++ b/matlab/gsa/mc_moments.m @@ -17,19 +17,23 @@ function [vdec, cc, ac] = mc_moments(mm, ss, dr) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global options_ M_ +global options_ M_ estim_params_ oo_ [nr1, nc1, nsam] = size(mm); + nobs=size(options_.varobs,1); disp('Computing theoretical moments ...') h = dyn_waitbar(0,'Theoretical moments ...'); + vdec = zeros(nobs,M_.exo_nbr,nsam); + cc = zeros(size(options_.varobs,1),size(options_.varobs,1),nsam); + ac = zeros(size(options_.varobs,1),size(options_.varobs,1)*options_.ar,nsam); for j=1:nsam, - dr.ghx = mm(:, [1:(nc1-M_.exo_nbr)],j); - dr.ghu = mm(:, [(nc1-M_.exo_nbr+1):end], j); + oo_.dr.ghx = mm(:, [1:(nc1-M_.exo_nbr)],j); + oo_.dr.ghu = mm(:, [(nc1-M_.exo_nbr+1):end], j); if ~isempty(ss), set_shocks_param(ss(j,:)); end - [vdec(:,:,j), corr, autocorr, z, zz] = th_moments(dr,options_.varobs); + [vdec(:,:,j), corr, autocorr, z, zz] = th_moments(oo_.dr,options_.varobs); cc(:,:,j)=triu(corr); dum=[]; for i=1:options_.ar diff --git a/matlab/gsa/th_moments.m b/matlab/gsa/th_moments.m index 2837fb696..9dbe7f8f6 100644 --- a/matlab/gsa/th_moments.m +++ b/matlab/gsa/th_moments.m @@ -35,11 +35,12 @@ function [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list) end end - [gamma_y,ivar] = th_autocovariances(dr,ivar,M_, options_); - m = dr.ys(ivar); + [gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_, options_); + m = dr.ys(ivar(stationary_vars)); - i1 = find(abs(diag(gamma_y{1})) > 1e-12); +% i1 = find(abs(diag(gamma_y{1})) > 1e-12); + i1 = [1:length(ivar)]; s2 = diag(gamma_y{1}); sd = sqrt(s2); From eb70ddfe66e4973e80a6f6062a4af33aae152932 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Tue, 10 Dec 2013 17:31:26 +0100 Subject: [PATCH 5/6] threshold_redform must override detailed mapping also with log-transformation --- matlab/gsa/redform_map.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/matlab/gsa/redform_map.m b/matlab/gsa/redform_map.m index 94c56a41d..2a2e9b027 100644 --- a/matlab/gsa/redform_map.m +++ b/matlab/gsa/redform_map.m @@ -132,7 +132,7 @@ for j=1:size(anamendo,1) iplo=iplo+1; js=js+1; xdir0 = [adir,filesep,namendo,'_vs_', namexo]; - if ilog==0, + if ilog==0 || ~isempty(threshold), if isempty(threshold) if isempty(dir(xdir0)) mkdir(xdir0) @@ -233,7 +233,7 @@ for j=1:size(anamendo,1) iplo=iplo+1; js=js+1; xdir0 = [adir,filesep,namendo,'_vs_', namlagendo]; - if ilog==0, + if ilog==0 || ~isempty(threshold), if isempty(threshold) if isempty(dir(xdir0)) mkdir(xdir0) From a90807f85bfab88c37db4d60d713ce3dec550029 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Tue, 10 Dec 2013 17:35:55 +0100 Subject: [PATCH 6/6] Cosmetic change --- matlab/gsa/mc_moments.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/matlab/gsa/mc_moments.m b/matlab/gsa/mc_moments.m index 4b0624145..5c2cd034a 100644 --- a/matlab/gsa/mc_moments.m +++ b/matlab/gsa/mc_moments.m @@ -24,8 +24,8 @@ global options_ M_ estim_params_ oo_ disp('Computing theoretical moments ...') h = dyn_waitbar(0,'Theoretical moments ...'); vdec = zeros(nobs,M_.exo_nbr,nsam); - cc = zeros(size(options_.varobs,1),size(options_.varobs,1),nsam); - ac = zeros(size(options_.varobs,1),size(options_.varobs,1)*options_.ar,nsam); + cc = zeros(nobs,nobs,nsam); + ac = zeros(nobs,nobs*options_.ar,nsam); for j=1:nsam, oo_.dr.ghx = mm(:, [1:(nc1-M_.exo_nbr)],j);