Merge pull request #561 from rattoma/master

Fixed bugs and continue GSA plots harmonization
time-shift
Sébastien Villemot 2013-12-11 01:51:50 -08:00
commit 1dc3f2119a
7 changed files with 79 additions and 71 deletions

View File

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

View File

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

View File

@ -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 <http://www.gnu.org/licenses/>.
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(nobs,nobs,nsam);
ac = zeros(nobs,nobs*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

View File

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

View File

@ -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_restrictions(Tt,Rr,M_,options_,oo_);
end
infox(j,1)=info(1);
if infox(j,1)==0 && ~exist('T'),
dr_=oo_.dr;
@ -468,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)';
@ -480,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_*.*']);
@ -533,6 +536,9 @@ if length(iunstable)>0 && length(iunstable)<Nsam,
if any(infox==30),
disp([' For ',num2str(length(find(infox==30))/Nsam*100,'%1.3f'),'\% Ergodic variance can''t be computed.'])
end
if any(infox==49),
disp([' For ',num2str(length(find(infox==49))/Nsam*100,'%1.3f'),'\% The model violates one (many) endogenous prior restriction(s).'])
end
end
skipline()
@ -597,7 +603,9 @@ if length(iunstable)>0 && length(iunstable)<Nsam,
c0=corrcoef(lpmat(istable,:));
c00=tril(c0,-1);
stab_map_2(lpmat(istable,:),alpha2, pvalue_corr, asname, OutputDirectoryName,xparam1,astitle);
if length(istable)>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

View File

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

View File

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