diff --git a/matlab/gsa/filt_mc_.m b/matlab/gsa/filt_mc_.m index 6cc14ce896d79c007ec753fe5fb93daaf3275905..2009299c38b29ee598a43fc14bbf37c611609d32 100644 --- a/matlab/gsa/filt_mc_.m +++ b/matlab/gsa/filt_mc_.m @@ -286,6 +286,7 @@ else PP(j,i)=P; end end + if ~options_.nograph, ifig=0; for i=1:size(vvarvecm,1), if mod(i,9)==1, @@ -376,6 +377,7 @@ else end end end + end param_names=''; for j=1:npar+nshock, @@ -470,6 +472,7 @@ else pnam=bayestopt_.name; % plot trade-offs + if ~options_.nograph a00=jet(size(vvarvecm,1)); for ix=1:ceil(length(nsnam)/5), hh = dyn_figure(options_); @@ -522,6 +525,7 @@ else end end close all + end for j=1:size(SP,2), nsx(j)=length(find(SP(:,j))); diff --git a/matlab/gsa/redform_map.m b/matlab/gsa/redform_map.m index 8fed92d98f312bc84792d51f86628607c3249d4f..f947432316c23dcb9d8e404b4be228eccaa8b492 100644 --- a/matlab/gsa/redform_map.m +++ b/matlab/gsa/redform_map.m @@ -14,7 +14,7 @@ function redform_map(dirname,options_gsa_) % Written by Marco Ratto % Joint Research Centre, The European Commission, % (http://eemc.jrc.ec.europa.eu/), -% marco.ratto@jrc.it +% marco.ratto@jrc.it % % Reference: % M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. @@ -56,26 +56,26 @@ pvalue_corr = options_gsa_.alpha2_redform; pnames = M_.param_names(estim_params_.param_vals(:,1),:); if nargin==0, - dirname=''; + dirname=''; end if pprior - load([dirname,'/',M_.fname,'_prior'],'lpmat', 'lpmat0', 'istable','T'); - adir=[dirname '/redform_stab']; + load([dirname,'/',M_.fname,'_prior'],'lpmat', 'lpmat0', 'istable','T'); + adir=[dirname '/redform_stab']; else - load([dirname,'/',M_.fname,'_mc'],'lpmat', 'lpmat0', 'istable','T'); - adir=[dirname '/redform_mc']; + load([dirname,'/',M_.fname,'_mc'],'lpmat', 'lpmat0', 'istable','T'); + adir=[dirname '/redform_mc']; end if ~exist('T') - stab_map_(dirname,options_gsa_); -if pprior - load([dirname,'/',M_.fname,'_prior'],'T'); -else - load([dirname,'/',M_.fname,'_mc'],'T'); -end + stab_map_(dirname,options_gsa_); + if pprior + load([dirname,'/',M_.fname,'_prior'],'T'); + else + load([dirname,'/',M_.fname,'_mc'],'T'); + end end if isempty(dir(adir)) - mkdir(adir) + mkdir(adir) end adir0=pwd; %cd(adir) @@ -86,242 +86,255 @@ xx0=lpmat0(istable,:); [kn, np]=size(x0); offset = length(bayestopt_.pshape)-np; if options_gsa_.prior_range, - pshape=5*(ones(np,1)); - pd = [NaN(np,1) NaN(np,1) bayestopt_.lb(offset+1:end) bayestopt_.ub(offset+1:end)]; + pshape=5*(ones(np,1)); + pd = [NaN(np,1) NaN(np,1) bayestopt_.lb(offset+1:end) bayestopt_.ub(offset+1:end)]; else - pshape = bayestopt_.pshape(offset+1:end); - pd = [bayestopt_.p6(offset+1:end) bayestopt_.p7(offset+1:end) bayestopt_.p3(offset+1:end) bayestopt_.p4(offset+1:end)]; + pshape = bayestopt_.pshape(offset+1:end); + pd = [bayestopt_.p6(offset+1:end) bayestopt_.p7(offset+1:end) bayestopt_.p3(offset+1:end) bayestopt_.p4(offset+1:end)]; end nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:))); -clear lpmat lpmat0 +clear lpmat lpmat0 js=0; for j=1:size(anamendo,1) - namendo=deblank(anamendo(j,:)); - iendo=strmatch(namendo,M_.endo_names(oo_.dr.order_var,:),'exact'); - ifig=0; - iplo=0; - for jx=1:size(anamexo,1) - namexo=deblank(anamexo(jx,:)); - iexo=strmatch(namexo,M_.exo_names,'exact'); + namendo=deblank(anamendo(j,:)); + iendo=strmatch(namendo,M_.endo_names(oo_.dr.order_var,:),'exact'); + ifig=0; + iplo=0; + for jx=1:size(anamexo,1) + namexo=deblank(anamexo(jx,:)); + iexo=strmatch(namexo,M_.exo_names,'exact'); + disp(' ') + disp(['[', namendo,' vs. ',namexo,']']) - if ~isempty(iexo), - %y0=squeeze(T(iendo,iexo+nspred,istable)); - y0=squeeze(T(iendo,iexo+nspred,:)); - if (max(y0)-min(y0))>1.e-10, - if mod(iplo,9)==0 && isempty(threshold), - ifig=ifig+1; - hfig = dyn_figure(options_,'name',[namendo,' vs. shocks ',int2str(ifig)]); - iplo=0; - end - iplo=iplo+1; - js=js+1; - xdir0 = [adir,'/',namendo,'_vs_', namexo]; - if ilog==0, - if isempty(threshold) - if isempty(dir(xdir0)) - mkdir(xdir0) - end - si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namexo, xdir0, options_gsa_); - else - iy=find( (y0>threshold(1)) & (y0<threshold(2))); - iyc=find( (y0<=threshold(1)) | (y0>=threshold(2))); - xdir = [xdir0,'_cut']; - if isempty(dir(xdir)) - mkdir(xdir) - end -% if ~isempty(iy), -% si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namexo, xdir, options_gsa_); -% else - si(:,js) = NaN(np,1); -% end - if ~isempty(iy) && ~isempty(iyc) - delete([xdir, '/*cut*.*']) - [proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0); -% indsmirnov = find(dproba>ksstat); - indsmirnov = find(proba<pvalue_ks); - stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir); - stab_map_2(x0(iy,:),alpha2,pvalue_corr,'cut',xdir) - stab_map_2(x0(iyc,:),alpha2,pvalue_corr,'trim',xdir) - lpmat=x0(iy,:); - lpmat0=xx0(iy,:); - istable=[1:length(iy)]; - save([xdir,filesep,'threshold.mat'],'lpmat','lpmat0','istable') - clear lpmat lpmat0 istable + + if ~isempty(iexo), + %y0=squeeze(T(iendo,iexo+nspred,istable)); + y0=squeeze(T(iendo,iexo+nspred,:)); + if (max(y0)-min(y0))>1.e-10, + if mod(iplo,9)==0 && isempty(threshold) && ~options_.nograph, + ifig=ifig+1; + hfig = dyn_figure(options_,'name',[namendo,' vs. shocks ',int2str(ifig)]); + iplo=0; + end + iplo=iplo+1; + js=js+1; + xdir0 = [adir,'/',namendo,'_vs_', namexo]; + if ilog==0, + if isempty(threshold) + if isempty(dir(xdir0)) + mkdir(xdir0) + end + si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namexo, xdir0, options_gsa_); + else + iy=find( (y0>threshold(1)) & (y0<threshold(2))); + iyc=find( (y0<=threshold(1)) | (y0>=threshold(2))); + xdir = [xdir0,'_cut']; + if isempty(dir(xdir)) + mkdir(xdir) + end + % if ~isempty(iy), + % si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namexo, xdir, options_gsa_); + % else + si(:,js) = NaN(np,1); + % end + if ~isempty(iy) && ~isempty(iyc) + delete([xdir, '/*cut*.*']) + [proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0); + % indsmirnov = find(dproba>ksstat); + indsmirnov = find(proba<pvalue_ks); + for jp=1:length(indsmirnov), + disp([M_.param_names(estim_params_.param_vals(indsmirnov(jp),1),:),' d-stat = ', num2str(dproba(indsmirnov(jp)),'%1.3f'),' p-value = ', num2str(proba(indsmirnov(jp)),'%1.3f')]) + end + disp(' '); + stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir); + stab_map_2(x0(iy,:),alpha2,pvalue_corr,'cut',xdir) + stab_map_2(x0(iyc,:),alpha2,pvalue_corr,'trim',xdir) + lpmat=x0(iy,:); + lpmat0=xx0(iy,:); + istable=[1:length(iy)]; + save([xdir,filesep,'threshold.mat'],'lpmat','lpmat0','istable') + clear lpmat lpmat0 istable + end + end + else + [yy, xdir] = log_trans_(y0,xdir0); + if isempty(dir(xdir)) + mkdir(xdir) + end + silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namexo, xdir, options_gsa_); + end + + if isempty(threshold) && ~options_.nograph, + figure(hfig) + subplot(3,3,iplo), + if ilog, + [saso, iso] = sort(-silog(:,js)); + bar([silog(iso(1:min(np,10)),js)]) + logflag='log'; + else + [saso, iso] = sort(-si(:,js)); + bar(si(iso(1:min(np,10)),js)) + logflag=''; + end + %set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8) + set(gca,'xticklabel',' ','fontsize',10) + set(gca,'xlim',[0.5 10.5]) + for ip=1:min(np,10), + text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') + end + title([logflag,' ',namendo,' vs. ',namexo],'interpreter','none') + if iplo==9, + dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_); + if ~options_.nodisplay + close(hfig); + end + end + end + end - end - else - [yy, xdir] = log_trans_(y0,xdir0); - if isempty(dir(xdir)) - mkdir(xdir) - end - silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namexo, xdir, options_gsa_); end - - if isempty(threshold) - figure(hfig) - subplot(3,3,iplo), - if ilog, - [saso, iso] = sort(-silog(:,js)); - bar([silog(iso(1:min(np,10)),js)]) - logflag='log'; - else - [saso, iso] = sort(-si(:,js)); - bar(si(iso(1:min(np,10)),js)) - logflag=''; - end - %set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8) - set(gca,'xticklabel',' ','fontsize',10) - set(gca,'xlim',[0.5 10.5]) - for ip=1:min(np,10), - text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') - end - title([logflag,' ',namendo,' vs. ',namexo],'interpreter','none') - if iplo==9, - dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_); - if ~options_.nodisplay + end + if iplo<9 && iplo>0 && ifig && ~options_.nograph, + dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_); + if ~options_.nodisplay close(hfig); - end end - end - - end - end - end - if iplo<9 && iplo>0 && ifig, - dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_); - if ~options_.nodisplay - close(hfig); end - end - ifig=0; - iplo=0; - for je=1:size(anamlagendo,1) - namlagendo=deblank(anamlagendo(je,:)); - ilagendo=strmatch(namlagendo,M_.endo_names(oo_.dr.order_var(oo_.dr.nstatic+1:oo_.dr.nstatic+nsok),:),'exact'); - - if ~isempty(ilagendo), - %y0=squeeze(T(iendo,ilagendo,istable)); - y0=squeeze(T(iendo,ilagendo,:)); - if (max(y0)-min(y0))>1.e-10, - if mod(iplo,9)==0 && isempty(threshold), - ifig=ifig+1; - hfig = dyn_figure(options_,'name',[namendo,' vs. lags ',int2str(ifig)]); - iplo=0; - end - iplo=iplo+1; - js=js+1; - xdir0 = [adir,'/',namendo,'_vs_', namlagendo]; - if ilog==0, - if isempty(threshold) - if isempty(dir(xdir0)) - mkdir(xdir0) - end - si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namlagendo, xdir0, options_gsa_); - else - iy=find( (y0>threshold(1)) & (y0<threshold(2))); - iyc=find( (y0<=threshold(1)) | (y0>=threshold(2))); - xdir = [xdir0,'_cut']; - if isempty(dir(xdir)) - mkdir(xdir) - end -% if ~isempty(iy) -% si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namlagendo, xdir, options_gsa_); -% end - if ~isempty(iy) && ~isempty(iyc), - delete([xdir, '/*cut*.*']) - [proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0); -% indsmirnov = find(dproba>ksstat); - indsmirnov = find(proba<pvalue_ks); - stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir); - stab_map_2(x0(iy,:),alpha2,pvalue_corr,'cut',xdir) - stab_map_2(x0(iyc,:),alpha2,pvalue_corr,'trim',xdir) - lpmat=x0(iy,:); - lpmat0=xx0(iy,:); - istable=[1:length(iy)]; - save([xdir,filesep,'threshold.mat'],'lpmat','lpmat0','istable') - clear lpmat lpmat0 istable - - end - end - else - [yy, xdir] = log_trans_(y0,xdir0); - if isempty(dir(xdir)) - mkdir(xdir) - end - silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namlagendo, xdir, options_gsa_); - end - - if isempty(threshold) - figure(hfig), - subplot(3,3,iplo), - if ilog, - [saso, iso] = sort(-silog(:,js)); - bar([silog(iso(1:min(np,10)),js)]) - logflag='log'; - else - [saso, iso] = sort(-si(:,js)); - bar(si(iso(1:min(np,10)),js)) - logflag=''; - end - %set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8) - set(gca,'xticklabel',' ','fontsize',10) - set(gca,'xlim',[0.5 10.5]) - for ip=1:min(np,10), - text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') - end - title([logflag,' ',namendo,' vs. ',namlagendo,'(-1)'],'interpreter','none') - if iplo==9, - dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_); - if ~options_.nodisplay - close(hfig); + ifig=0; + iplo=0; + for je=1:size(anamlagendo,1) + namlagendo=deblank(anamlagendo(je,:)); + ilagendo=strmatch(namlagendo,M_.endo_names(oo_.dr.order_var(oo_.dr.nstatic+1:oo_.dr.nstatic+nsok),:),'exact'); + disp(' ') + disp(['[', namendo,' vs. lagged ',namlagendo,']']) + + if ~isempty(ilagendo), + %y0=squeeze(T(iendo,ilagendo,istable)); + y0=squeeze(T(iendo,ilagendo,:)); + if (max(y0)-min(y0))>1.e-10, + if mod(iplo,9)==0 && isempty(threshold) && ~options_.nograph, + ifig=ifig+1; + hfig = dyn_figure(options_,'name',[namendo,' vs. lags ',int2str(ifig)]); + iplo=0; + end + iplo=iplo+1; + js=js+1; + xdir0 = [adir,'/',namendo,'_vs_', namlagendo]; + if ilog==0, + if isempty(threshold) + if isempty(dir(xdir0)) + mkdir(xdir0) + end + si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namlagendo, xdir0, options_gsa_); + else + iy=find( (y0>threshold(1)) & (y0<threshold(2))); + iyc=find( (y0<=threshold(1)) | (y0>=threshold(2))); + xdir = [xdir0,'_cut']; + if isempty(dir(xdir)) + mkdir(xdir) + end + % if ~isempty(iy) + % si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namlagendo, xdir, options_gsa_); + % end + if ~isempty(iy) && ~isempty(iyc), + delete([xdir, '/*cut*.*']) + [proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0); + % indsmirnov = find(dproba>ksstat); + indsmirnov = find(proba<pvalue_ks); + for jp=1:length(indsmirnov), + disp([M_.param_names(estim_params_.param_vals(indsmirnov(jp),1),:),' d-stat = ', num2str(dproba(indsmirnov(jp)),'%1.3f'),' p-value = ', num2str(proba(indsmirnov(jp)),'%1.3f')]) + end + disp(' '); + stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir); + stab_map_2(x0(iy,:),alpha2,pvalue_corr,'cut',xdir) + stab_map_2(x0(iyc,:),alpha2,pvalue_corr,'trim',xdir) + lpmat=x0(iy,:); + lpmat0=xx0(iy,:); + istable=[1:length(iy)]; + save([xdir,filesep,'threshold.mat'],'lpmat','lpmat0','istable') + clear lpmat lpmat0 istable + + end + end + else + [yy, xdir] = log_trans_(y0,xdir0); + if isempty(dir(xdir)) + mkdir(xdir) + end + silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namlagendo, xdir, options_gsa_); + end + + if isempty(threshold) && ~options_.nograph + figure(hfig), + subplot(3,3,iplo), + if ilog, + [saso, iso] = sort(-silog(:,js)); + bar([silog(iso(1:min(np,10)),js)]) + logflag='log'; + else + [saso, iso] = sort(-si(:,js)); + bar(si(iso(1:min(np,10)),js)) + logflag=''; + end + %set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8) + set(gca,'xticklabel',' ','fontsize',10) + set(gca,'xlim',[0.5 10.5]) + for ip=1:min(np,10), + text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') + end + title([logflag,' ',namendo,' vs. ',namlagendo,'(-1)'],'interpreter','none') + if iplo==9, + dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_); + if ~options_.nodisplay + close(hfig); + end + end + end + end end - end - - end end - end - if iplo<9 && iplo>0 && ifig, - dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_); - if ~options_.nodisplay - close(hfig); + if iplo<9 && iplo>0 && ifig && ~options_.nograph, + dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_); + if ~options_.nodisplay + close(hfig); + end end - end end -if isempty(threshold), -if ilog==0, -hfig=dyn_figure(options_); %bar(si) -% boxplot(si','whis',10,'symbol','r.') -myboxplot(si',[],'.',[],10) -xlabel(' ') -set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np]) -set(gca,'xlim',[0.5 np+0.5]) -set(gca,'ylim',[0 1]) -set(gca,'position',[0.13 0.2 0.775 0.7]) -for ip=1:np, - text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') -end -title('Reduced form GSA') -dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_gsa'],options_); - -else -hfig=dyn_figure(options_); %bar(silog) -% boxplot(silog','whis',10,'symbol','r.') -myboxplot(silog',[],'.',[],10) -set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np]) -xlabel(' ') -set(gca,'xlim',[0.5 np+0.5]) -set(gca,'ylim',[0 1]) -set(gca,'position',[0.13 0.2 0.775 0.7]) -for ip=1:np, - text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') -end -title('Reduced form GSA - Log-transformed elements') -dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_gsa_log'],options_); - -end +if isempty(threshold) && ~options_.nograph, + if ilog==0, + hfig=dyn_figure(options_); %bar(si) + % boxplot(si','whis',10,'symbol','r.') + myboxplot(si',[],'.',[],10) + xlabel(' ') + set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np]) + set(gca,'xlim',[0.5 np+0.5]) + set(gca,'ylim',[0 1]) + set(gca,'position',[0.13 0.2 0.775 0.7]) + for ip=1:np, + text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') + end + title('Reduced form GSA') + dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_gsa'],options_); + + else + hfig=dyn_figure(options_); %bar(silog) + % boxplot(silog','whis',10,'symbol','r.') + myboxplot(silog',[],'.',[],10) + set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np]) + xlabel(' ') + set(gca,'xlim',[0.5 np+0.5]) + set(gca,'ylim',[0 1]) + set(gca,'position',[0.13 0.2 0.775 0.7]) + for ip=1:np, + text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') + end + title('Reduced form GSA - Log-transformed elements') + dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_gsa_log'],options_); + + end end function si = redform_private(x0, y0, pshape, pd, iload, pnames, namy, namx, xdir, opt_gsa) @@ -330,69 +343,78 @@ global bayestopt_ options_ % opt_gsa=options_.opt_gsa; np=size(x0,2); x00=x0; - if opt_gsa.prior_range, +if opt_gsa.prior_range, for j=1:np, - x0(:,j)=(x0(:,j)-pd(j,3))./(pd(j,4)-pd(j,3)); + x0(:,j)=(x0(:,j)-pd(j,3))./(pd(j,4)-pd(j,3)); end - else +else x0=priorcdf(x0,pshape, pd(:,1), pd(:,2), pd(:,3), pd(:,4)); - end +end fname=[xdir,'/map']; if iload==0, - hfig=dyn_figure(options_); hist(y0,30), title([namy,' vs. ', namx]) - if isempty(dir(xdir)) - mkdir(xdir) - end - dyn_saveas(hfig,[xdir,'/', namy,'_vs_', namx],options_); - if ~options_.nodisplay - close(hfig); - end -% gsa_ = gsa_sdp_dyn(y0, x0, -2, [],[],[],1,fname, pnames); - nrun=length(y0); - nest=min(250,nrun); - nfit=min(1000,nrun); -% dotheplots = (nfit<=nest); - gsa_ = gsa_sdp(y0(1:nest), x0(1:nest,:), 2, [],[-1 -1 -1 -1 -1 0],[],0,[fname,'_est'], pnames); - if nfit>nest, - gsa_ = gsa_sdp(y0(1:nfit), x0(1:nfit,:), -2, gsa_.nvr*nest^3/nfit^3,[-1 -1 -1 -1 -1 0],[],0,fname, pnames); - end - save([fname,'.mat'],'gsa_') - [sidum, iii]=sort(-gsa_.si); - gsa_.x0=x00(1:nfit,:); - hfig=gsa_sdp_plot(gsa_,fname,pnames,iii(1:min(12,np))); - close(hfig); - gsa_.x0=x0(1:nfit,:); -% copyfile([fname,'_est.mat'],[fname,'.mat']) - hfig=dyn_figure(options_); - plot(y0(1:nfit),[gsa_.fit y0(1:nfit)],'.'), - title([namy,' vs. ', namx,' fit']) - dyn_saveas(hfig,[xdir,'/', namy,'_vs_', namx,'_fit'],options_); - if ~options_.nodisplay - close(hfig); - end - if nfit<nrun, - npred=[nfit+1:nrun]; - yf = ss_anova_fcast(x0(npred,:), gsa_); - hfig=dyn_figure(options_); - plot(y0(npred),[yf y0(npred)],'.'), - title([namy,' vs. ', namx,' pred']) - dyn_saveas(hfig,[xdir,'/', namy,'_vs_', namx,'_pred'],options_); - if ~options_.nodisplay - close(hfig); - end - end + if isempty(dir(xdir)) + mkdir(xdir) + end + if ~options_.nograph, + hfig=dyn_figure(options_); hist(y0,30), title([namy,' vs. ', namx]) + dyn_saveas(hfig,[xdir,'/', namy,'_vs_', namx],options_); + if ~options_.nodisplay + close(hfig); + end + end + % gsa_ = gsa_sdp_dyn(y0, x0, -2, [],[],[],1,fname, pnames); + nrun=length(y0); + nest=min(250,nrun); + nfit=min(1000,nrun); + % dotheplots = (nfit<=nest); + gsa_ = gsa_sdp(y0(1:nest), x0(1:nest,:), 2, [],[-1 -1 -1 -1 -1 0],[],0,[fname,'_est'], pnames); + if nfit>nest, + gsa_ = gsa_sdp(y0(1:nfit), x0(1:nfit,:), -2, gsa_.nvr*nest^3/nfit^3,[-1 -1 -1 -1 -1 0],[],0,fname, pnames); + end + save([fname,'.mat'],'gsa_') + [sidum, iii]=sort(-gsa_.si); + gsa_.x0=x00(1:nfit,:); + if ~options_.nograph, + hfig=gsa_sdp_plot(gsa_,fname,pnames,iii(1:min(12,np))); + close(hfig); + end + gsa_.x0=x0(1:nfit,:); + % copyfile([fname,'_est.mat'],[fname,'.mat']) + if ~options_.nograph, + hfig=dyn_figure(options_); + plot(y0(1:nfit),[gsa_.fit y0(1:nfit)],'.'), + title([namy,' vs. ', namx,' fit']) + dyn_saveas(hfig,[xdir,'/', namy,'_vs_', namx,'_fit'],options_); + if ~options_.nodisplay + close(hfig); + end + if nfit<nrun, + npred=[nfit+1:nrun]; + yf = ss_anova_fcast(x0(npred,:), gsa_); + hfig=dyn_figure(options_); + plot(y0(npred),[yf y0(npred)],'.'), + title([namy,' vs. ', namx,' pred']) + dyn_saveas(hfig,[xdir,'/', namy,'_vs_', namx,'_pred'],options_); + if ~options_.nodisplay + close(hfig); + end + end + + end else -% gsa_ = gsa_sdp_dyn(y0, x0, 0, [],[],[],0,fname, pnames); - gsa_ = gsa_sdp(y0, x0, 0, [],[],[],0,fname, pnames); - yf = ss_anova_fcast(x0, gsa_); - hfig=dyn_figure(options_); - plot(y0,[yf y0],'.'), - title([namy,' vs. ', namx,' pred']) - dyn_saveas(hfig,[xdir,'/', namy,'_vs_', namx,'_pred'],options_); - if ~options_.nodisplay - close(hfig); - end + % gsa_ = gsa_sdp_dyn(y0, x0, 0, [],[],[],0,fname, pnames); + gsa_ = gsa_sdp(y0, x0, 0, [],[],[],0,fname, pnames); + if ~options_.nograph, + yf = ss_anova_fcast(x0, gsa_); + hfig=dyn_figure(options_); + plot(y0,[yf y0],'.'), + title([namy,' vs. ', namx,' pred']) + dyn_saveas(hfig,[xdir,'/', namy,'_vs_', namx,'_pred'],options_); + if ~options_.nodisplay + close(hfig); + end + end end % si = gsa_.multivariate.si; si = gsa_.si; diff --git a/matlab/gsa/stab_map_1.m b/matlab/gsa/stab_map_1.m index 3c354cc994f6469266e0040ee44290e830ee1f10..85df2d838b603da9c5e6c4786453eab6418adf90 100644 --- a/matlab/gsa/stab_map_1.m +++ b/matlab/gsa/stab_map_1.m @@ -77,7 +77,7 @@ if isempty(ipar), ipar=find(proba<pcrit); end nparplot=length(ipar); -if iplot +if iplot && ~options_.nograph lpmat=lpmat(:,ipar); ftit=bayestopt_.name(ipar+nshock*(1-ishock)); diff --git a/matlab/gsa/stab_map_2.m b/matlab/gsa/stab_map_2.m index aecfabf81e00c5daea1dd25b85f70a770e5c5426..b4de5b8cb3b445f218a87fb048d55848b2088c42 100644 --- a/matlab/gsa/stab_map_2.m +++ b/matlab/gsa/stab_map_2.m @@ -64,57 +64,70 @@ if ishock==0 else npar=estim_params_.np+nshock; end +disp([' ']) +disp(['Correlation analysis for ',fnam]) + for j=1:npar, - i2=find(abs(c00(:,j))>alpha2); - if length(i2)>0, - for jx=1:length(i2), - tval = abs(c00(i2(jx),j)*sqrt( (nsam-2)/(1-c00(i2(jx),j)^2) )); - tcr = tcrit(nsam-2,pvalue); - if tval>tcr, - j2=j2+1; - if mod(j2,12)==1, - ifig=ifig+1; - hh=dyn_figure(options_,'name',['Correlations in the ',fnam,' sample ', num2str(ifig)]); - end - subplot(3,4,j2-(ifig-1)*12) - % bar(c0(i2,j)), - % set(gca,'xticklabel',bayestopt_.name(i2)), - % set(gca,'xtick',[1:length(i2)]) - %plot(stock_par(ixx(nfilt+1:end,i),j),stock_par(ixx(nfilt+1:end,i),i2(jx)),'.k') - %hold on, - plot(x(:,j),x(:,i2(jx)),'.') - if ~isempty(xparam1) - hold on, plot(xparam1(j),xparam1(i2(jx)),'ro') - end - % xlabel(deblank(estim_params_.param_names(j,:)),'interpreter','none'), - % ylabel(deblank(estim_params_.param_names(i2(jx),:)),'interpreter','none'), - if ishock, - xlabel(bayestopt_.name{j},'interpreter','none'), - ylabel(bayestopt_.name{i2(jx)},'interpreter','none'), - else - xlabel(bayestopt_.name{j+nshock},'interpreter','none'), - ylabel(bayestopt_.name{i2(jx)+nshock},'interpreter','none'), - end - title(['cc = ',num2str(c0(i2(jx),j))]) - if (mod(j2,12)==0) && j2>0, + i2=find(abs(c00(:,j))>alpha2); + if length(i2)>0, + for jx=1:length(i2), + tval = abs(c00(i2(jx),j)*sqrt( (nsam-2)/(1-c00(i2(jx),j)^2) )); + tcr = tcrit(nsam-2,pvalue); + if tval>tcr, + j2=j2+1; + if ishock, + tmp_name = (['[',bayestopt_.name{j},',',bayestopt_.name{i2(jx)},']']); + else + tmp_name = (['[',bayestopt_.name{j+nshock},',',bayestopt_.name{i2(jx)+nshock},']']); + end + fprintf(1,'%20s: corrcoef = %7.3f\n',tmp_name,c0(i2(jx),j)); + + if ~options_.nograph, + + if mod(j2,12)==1, + ifig=ifig+1; + hh=dyn_figure(options_,'name',['Correlations in the ',fnam,' sample ', num2str(ifig)]); + end + subplot(3,4,j2-(ifig-1)*12) + % bar(c0(i2,j)), + % set(gca,'xticklabel',bayestopt_.name(i2)), + % set(gca,'xtick',[1:length(i2)]) + %plot(stock_par(ixx(nfilt+1:end,i),j),stock_par(ixx(nfilt+1:end,i),i2(jx)),'.k') + %hold on, + plot(x(:,j),x(:,i2(jx)),'.') + if ~isempty(xparam1) + hold on, plot(xparam1(j),xparam1(i2(jx)),'ro') + end + % xlabel(deblank(estim_params_.param_names(j,:)),'interpreter','none'), + % ylabel(deblank(estim_params_.param_names(i2(jx),:)),'interpreter','none'), + if ishock, + xlabel(bayestopt_.name{j},'interpreter','none'), + ylabel(bayestopt_.name{i2(jx)},'interpreter','none'), + else + xlabel(bayestopt_.name{j+nshock},'interpreter','none'), + ylabel(bayestopt_.name{i2(jx)+nshock},'interpreter','none'), + end + title(['cc = ',num2str(c0(i2(jx),j))]) + if (mod(j2,12)==0) && j2>0, + dyn_saveas(hh,[dirname,'/',fig_nam_,int2str(ifig)],options_); + if ~options_.nodisplay + close(hh); + end + end + end + end + + end + end + if ~options_.nograph && (j==(npar)) && j2>0, dyn_saveas(hh,[dirname,'/',fig_nam_,int2str(ifig)],options_); if ~options_.nodisplay - close(hh); + close(hh); end - end - end - end - end - if (j==(npar)) && j2>0, - dyn_saveas(hh,[dirname,'/',fig_nam_,int2str(ifig)],options_); - if ~options_.nodisplay - close(hh); - end - end - + end -if ifig==0, +if j2==0, disp(['No correlation term with pvalue <', num2str(pvalue),' found for ',fnam]) end %close all