Merge branch 'marco_gsa_fixes'

time-shift
Stéphane Adjemian (Charybdis) 2015-04-14 15:54:59 +02:00
commit da2c529ee7
7 changed files with 439 additions and 184 deletions

View File

@ -163,6 +163,11 @@ options_gsa = set_default_option(options_gsa,'istart_rmse',options_.presample+1)
options_gsa = set_default_option(options_gsa,'alpha_rmse',0.001);
options_gsa = set_default_option(options_gsa,'alpha2_rmse',1.e-5);
if options_gsa.neighborhood_width,
options_gsa.pprior=0;
options_gsa.ppost=0;
end
if options_gsa.redform && options_gsa.neighborhood_width==0 && isempty(options_gsa.threshold_redform),
options_gsa.pprior=1;
options_gsa.ppost=0;

View File

@ -43,7 +43,7 @@ alpha = options_gsa_.alpha_rmse;
% alpha2 = options_gsa_.alpha2_rmse;
alpha2 = 0;
pvalue = options_gsa_.alpha2_rmse;
istart = options_gsa_.istart_rmse;
istart = max(2,options_gsa_.istart_rmse);
alphaPC = 0.5;
fname_ = M_.fname;
@ -255,28 +255,32 @@ if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
stab_map_1(x, ilik(1:nfilt), ilik(nfilt+1:end), anam, 1,[],OutDir);
stab_map_2(x(ilik(1:nfilt),:),alpha2,pvalue,anam, OutDir);
else
for i=1:size(vvarvecm,1),
[dum, ixx(:,i)]=sort(rmse_MC(:,i));
if options_.opt_gsa.ppost,
if options_.opt_gsa.ppost,
rmse_txt=rmse_pmean;
r2_txt=r2_pmean;
else
if options_.opt_gsa.pprior || ~exist('rmse_pmean'),
if exist('rmse_mode'),
rmse_txt=rmse_mode;
r2_txt=r2_mode;
else
rmse_txt=NaN(1,size(rmse_MC,2));
r2_txt=NaN(1,size(r2_MC,2));
end
else
%nfilt0(i)=length(find(rmse_MC(:,i)<rmse_pmean(i)));
rmse_txt=rmse_pmean;
r2_txt=r2_pmean;
else
if options_.opt_gsa.pprior || ~exist('rmse_pmean'),
if exist('rmse_mode'),
rmse_txt=rmse_mode;
r2_txt=r2_mode;
else
rmse_txt=NaN(1,size(rmse_MC,2));
r2_txt=NaN(1,size(r2_MC,2));
end
else
%nfilt0(i)=length(find(rmse_MC(:,i)<rmse_pmean(i)));
rmse_txt=rmse_pmean;
r2_txt=r2_pmean;
end
end
for j=1:npar+nshock,
end
for i=1:size(vvarvecm,1),
[dum, ixx(:,i)]=sort(rmse_MC(:,i));
end
PP=ones(npar+nshock,size(vvarvecm,1));
PPV=ones(size(vvarvecm,1),size(vvarvecm,1),npar+nshock);
SS=zeros(npar+nshock,size(vvarvecm,1));
for j=1:npar+nshock,
for i=1:size(vvarvecm,1),
[H,P,KSSTAT] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j), alpha);
[H1,P1,KSSTAT1] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,1);
[H2,P2,KSSTAT2] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,-1);
@ -289,6 +293,22 @@ else
end
PP(j,i)=P;
end
for i=1:size(vvarvecm,1),
for l=1:size(vvarvecm,1),
if l~=i && PP(j,i)<alpha && PP(j,l)<alpha
[H,P,KSSTAT] = smirnov(x(ixx(1:nfilt0(i),i),j),x(ixx(1:nfilt0(l),l),j), alpha);
%[H1,P1,KSSTAT1] = smirnov(x(ixx(1:nfilt0(i),i),j),x(:,j), alpha);
% PP(j,i)=min(P,PP(j,i));
% PP(j,l)=min(P,PP(j,l));
%if P<P1
% if SS(j,i)*SS(j,l)
PPV(i,l,j) = P;
% end
elseif l==i
PPV(i,l,j) = PP(j,i);
end
end
end
end
if ~options_.nograph,
ifig=0;
@ -308,12 +328,17 @@ else
end
subplot(3,3,i-9*(ifig-1))
h=cumplot(lnprior(ixx(1:nfilt0(i),i)));
set(h,'color','red')
hold on, cumplot(lnprior)
set(h,'color','blue','linewidth',2)
hold on, h=cumplot(lnprior);
set(h,'color','k','linewidth',1)
h=cumplot(lnprior(ixx(nfilt0(i)+1:end,i)));
set(h,'color','green')
set(h,'color','red','linewidth',2)
title(vvarvecm(i,:),'interpreter','none')
if mod(i,9)==0 || i==size(vvarvecm,1)
if ~isoctave
annotation('textbox', [0.1,0,0.35,0.05],'String', 'Log-prior for BETTER R2','Color','Blue','horizontalalignment','center');
annotation('textbox', [0.55,0,0.35,0.05],'String', 'Log-prior for WORSE R2', 'Color','Red','horizontalalignment','center');
end
if options_.opt_gsa.ppost
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_post_lnprior',int2str(ifig)],options_);
else
@ -342,15 +367,20 @@ else
end
subplot(3,3,i-9*(ifig-1))
h=cumplot(likelihood(ixx(1:nfilt0(i),i)));
set(h,'color','red')
set(h,'color','blue','linewidth',2)
hold on, h=cumplot(likelihood);
set(h,'color','k','linewidth',1)
h=cumplot(likelihood(ixx(nfilt0(i)+1:end,i)));
set(h,'color','green')
set(h,'color','red','linewidth',2)
title(vvarvecm(i,:),'interpreter','none')
if options_.opt_gsa.ppost==0,
set(gca,'xlim',[min( likelihood(ixx(1:nfilt0(i),i)) ) max( likelihood(ixx(1:nfilt0(i),i)) )])
end
if mod(i,9)==0 || i==size(vvarvecm,1)
if ~isoctave
annotation('textbox', [0.1,0,0.35,0.05],'String', 'Log-likelihood for BETTER R2','Color','Blue','horizontalalignment','center');
annotation('textbox', [0.55,0,0.35,0.05],'String', 'Log-likelihood for WORSE R2', 'Color','Red','horizontalalignment','center');
end
if options_.opt_gsa.ppost
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_post_lnlik',int2str(ifig) ],options_);
else
@ -379,15 +409,20 @@ else
end
subplot(3,3,i-9*(ifig-1))
h=cumplot(logpo2(ixx(1:nfilt0(i),i)));
set(h,'color','red')
set(h,'color','blue','linewidth',2)
hold on, h=cumplot(logpo2);
set(h,'color','k','linewidth',1)
h=cumplot(logpo2(ixx(nfilt0(i)+1:end,i)));
set(h,'color','green')
set(h,'color','red','linewidth',2)
title(vvarvecm(i,:),'interpreter','none')
if options_.opt_gsa.ppost==0,
set(gca,'xlim',[min( logpo2(ixx(1:nfilt0(i),i)) ) max( logpo2(ixx(1:nfilt0(i),i)) )])
end
if mod(i,9)==0 || i==size(vvarvecm,1)
if ~isoctave
annotation('textbox', [0.1,0,0.35,0.05],'String', 'Log-posterior for BETTER R2','Color','Blue','horizontalalignment','center');
annotation('textbox', [0.55,0,0.35,0.05],'String', 'Log-posterior for WORSE R2', 'Color','Red','horizontalalignment','center');
end
if options_.opt_gsa.ppost
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_post_lnpost',int2str(ifig) ],options_);
else
@ -401,11 +436,12 @@ else
end
end
param_names='';
for j=1:npar+nshock,
param_names=char(param_names, bayestopt_.name{j});
end
param_names=param_names(2:end,:);
param_names=char(bayestopt_.name);
% param_names='';
% for j=1:npar+nshock,
% param_names=char(param_names, bayestopt_.name{j});
% end
% param_names=param_names(2:end,:);
skipline()
disp('RMSE over the MC sample:')
@ -526,15 +562,94 @@ else
a00=jet(size(vvarvecm,1));
if options_.opt_gsa.ppost
temp_name='RMSE Posterior Tradeoffs:';
atitle='RMSE Posterior Map:';
asname='rmse_post';
else
if options_.opt_gsa.pprior
temp_name='RMSE Prior Tradeoffs:';
atitle='RMSE Prior Map:';
asname='rmse_prior';
else
temp_name='RMSE MC Tradeoffs;';
temp_name='RMSE MC Tradeoffs:';
atitle='RMSE MC Map:';
asname='rmse_mc';
end
end
% now I plot by observed variables
options_mcf.pvalue_ks = alpha;
options_mcf.pvalue_corr = pvalue;
options_mcf.alpha2 = alpha2;
options_mcf.param_names = char(bayestopt_.name);
options_mcf.fname_ = fname_;
options_mcf.OutputDirectoryName = OutDir;
for iy=1:size(vvarvecm,1),
options_mcf.amcf_name = [asname '_' deblank(vvarvecm(iy,:)) '_map' ];
options_mcf.amcf_title = [atitle ' ' deblank(vvarvecm(iy,:))];
options_mcf.beha_title = ['better fit of ' deblank(vvarvecm(iy,:))];
options_mcf.nobeha_title = ['worse fit of ' deblank(vvarvecm(iy,:))];
options_mcf.title = ['the fit of ' deblank(vvarvecm(iy,:))];
mcf_analysis(x, ixx(1:nfilt0(iy),iy), ixx(nfilt0(iy)+1:end,iy), options_mcf, options_)
end
for iy=1:size(vvarvecm,1),
ipar = find(any(squeeze(PPV(iy,:,:))<alpha));
for ix=1:ceil(length(ipar)/5),
hh = dyn_figure(options_,'name',[temp_name,' observed variable ',deblank(vvarvecm(iy,:))]);
for j=1+5*(ix-1):min(length(ipar),5*ix),
subplot(2,3,j-5*(ix-1))
%h0=cumplot(x(:,nsnam(j)+nshock));
h0=cumplot(x(:,ipar(j)));
set(h0,'color',[0 0 0])
hold on,
iobs=find(squeeze(PPV(iy,:,ipar(j)))<alpha);
for i=1:size(vvarvecm,1),
%h0=cumplot(x(ixx(1:nfilt,np(i)),nsnam(j)+nshock));
% h0=cumplot(x(ixx(1:nfilt0(np(i)),np(i)),nsnam(j)));
if any(iobs==i) || i==iy,
h0=cumplot(x(ixx(1:nfilt0(i),i),ipar(j)));
if ~isoctave
hcmenu = uicontextmenu;
uimenu(hcmenu,'Label',deblank(vvarvecm(i,:)));
set(h0,'uicontextmenu',hcmenu)
end
else
h0=cumplot(x(ixx(1:nfilt0(i),i),ipar(j))*NaN);
end
set(h0,'color',a00(i,:),'linewidth',2)
end
ydum=get(gca,'ylim');
%xdum=xparam1(nshock+nsnam(j));
if exist('xparam1')
xdum=xparam1(ipar(j));
h1=plot([xdum xdum],ydum);
set(h1,'color',[0.85 0.85 0.85],'linewidth',2)
end
xlabel('')
title([pnam{ipar(j)}],'interpreter','none')
end
%subplot(3,2,6)
if isoctave
legend(char('base',vvarvecm),'location','eastoutside');
else
h0=legend(char('base',vvarvecm),0);
set(h0,'fontsize',6,'position',[0.7 0.1 0.2 0.3],'interpreter','none');
end
%h0=legend({'base',vnam{np}}',0);
%set(findobj(get(h0,'children'),'type','text'),'interpreter','none')
if options_.opt_gsa.ppost
dyn_saveas(hh,[ OutDir filesep fname_ '_rmse_post_' deblank(vvarvecm(iy,:)) '_' int2str(ix)],options_);
else
if options_.opt_gsa.pprior
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_prior_' deblank(vvarvecm(iy,:)) '_' int2str(ix) ],options_);
else
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_mc_' deblank(vvarvecm(iy,:)) '_' int2str(ix)],options_);
end
end
end
end
% now I plot by individual parameters
for ix=1:ceil(length(nsnam)/5),
hh = dyn_figure(options_,'name',[temp_name,' observed variables ',int2str(ix)]);
hh = dyn_figure(options_,'name',[temp_name,' estimated params and shocks ',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));
@ -551,8 +666,13 @@ else
h0=cumplot(x(ixx(1:nfilt0(i),i),nsnam(j))*NaN);
else
h0=cumplot(x(ixx(1:nfilt0(i),i),nsnam(j)));
end
set(h0,'color',a00(i,:))
if ~isoctave
hcmenu = uicontextmenu;
uimenu(hcmenu,'Label',deblank(vvarvecm(i,:)));
set(h0,'uicontextmenu',hcmenu)
end
end
set(h0,'color',a00(i,:),'linewidth',2)
end
ydum=get(gca,'ylim');
%xdum=xparam1(nshock+nsnam(j));
@ -574,74 +694,74 @@ else
%h0=legend({'base',vnam{np}}',0);
%set(findobj(get(h0,'children'),'type','text'),'interpreter','none')
if options_.opt_gsa.ppost
dyn_saveas(hh,[ OutDir filesep fname_ '_rmse_post_' int2str(ix)],options_);
dyn_saveas(hh,[ OutDir filesep fname_ '_rmse_post_params_' int2str(ix)],options_);
else
if options_.opt_gsa.pprior
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_prior_' int2str(ix) ],options_);
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_prior_params_' int2str(ix) ],options_);
else
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_mc_' int2str(ix)],options_);
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_mc_params_' int2str(ix)],options_);
end
end
end
end
for j=1:size(SP,2),
nsx(j)=length(find(SP(:,j)));
end
% for j=1:size(SP,2),
% nsx(j)=length(find(SP(:,j)));
% end
skipline(2)
disp('Sensitivity table (significance and direction):')
vav=char(zeros(1, size(param_names,2)+3 ));
ibl = 12-size(vvarvecm,2);
for j=1:size(vvarvecm,1),
vav = [vav, char(zeros(1,ibl)),vvarvecm(j,:)];
end
disp(vav)
for j=1:npar+nshock, %estim_params_.np,
%disp([param_names(j,:), sprintf('%8.5g',SP(j,:))])
disp([param_names(j,:),' ', sprintf('%12.3g',PP(j,:))])
disp([char(zeros(1, size(param_names,2)+3 )),sprintf(' (%6g)',SS(j,:))])
end
% skipline(2)
% disp('Sensitivity table (significance and direction):')
% vav=char(zeros(1, size(param_names,2)+3 ));
% ibl = 12-size(vvarvecm,2);
% for j=1:size(vvarvecm,1),
% vav = [vav, char(zeros(1,ibl)),vvarvecm(j,:)];
% end
% disp(vav)
% for j=1:npar+nshock, %estim_params_.np,
% %disp([param_names(j,:), sprintf('%8.5g',SP(j,:))])
% disp([param_names(j,:),' ', sprintf('%12.3g',PP(j,:))])
% disp([char(zeros(1, size(param_names,2)+3 )),sprintf(' (%6g)',SS(j,:))])
% end
skipline()
disp('Starting bivariate analysis:')
for i=1:size(vvarvecm,1)
if options_.opt_gsa.ppost
fnam = ['rmse_post_',deblank(vvarvecm(i,:))];
else
if options_.opt_gsa.pprior
fnam = ['rmse_prior_',deblank(vvarvecm(i,:))];
else
fnam = ['rmse_mc_',deblank(vvarvecm(i,:))];
end
end
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)])
% ifig=0;
% j2=0;
% for j=1:npar+nshock,
% i2=find(abs(pc(:,j))>alphaPC);
% if ~isempty(i2),
% j2=j2+1;
% if mod(j2,12)==1,
% ifig=ifig+1;
% figure('name',['PCA of the filtered sample ',deblank(vvarvecm(i,:)),' ',num2str(ifig)]),
% end
% subplot(3,4,j2-(ifig-1)*12)
% bar(pc(i2,j)),
% set(gca,'xticklabel',bayestopt_.name(i2)),
% set(gca,'xtick',[1:length(i2)])
% title(['PC ',num2str(j),'. Explained ',num2str(explained(j)),'%'])
% end
% if (mod(j2,12)==0 | j==(npar+nshock)) & j2,
% saveas(gcf,[fname_,'_SA_PCA_',deblank(vvarvecm(i,:)),'_',int2str(ifig)])
% end
% end
% close all
end
% skipline()
% disp('Starting bivariate analysis:')
%
% for i=1:size(vvarvecm,1)
% if options_.opt_gsa.ppost
% fnam = ['rmse_post_',deblank(vvarvecm(i,:))];
% else
% if options_.opt_gsa.pprior
% fnam = ['rmse_prior_',deblank(vvarvecm(i,:))];
% else
% fnam = ['rmse_mc_',deblank(vvarvecm(i,:))];
% end
% end
% 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)])
% % ifig=0;
% % j2=0;
% % for j=1:npar+nshock,
% % i2=find(abs(pc(:,j))>alphaPC);
% % if ~isempty(i2),
% % j2=j2+1;
% % if mod(j2,12)==1,
% % ifig=ifig+1;
% % figure('name',['PCA of the filtered sample ',deblank(vvarvecm(i,:)),' ',num2str(ifig)]),
% % end
% % subplot(3,4,j2-(ifig-1)*12)
% % bar(pc(i2,j)),
% % set(gca,'xticklabel',bayestopt_.name(i2)),
% % set(gca,'xtick',[1:length(i2)])
% % title(['PC ',num2str(j),'. Explained ',num2str(explained(j)),'%'])
% % end
% % if (mod(j2,12)==0 | j==(npar+nshock)) & j2,
% % saveas(gcf,[fname_,'_SA_PCA_',deblank(vvarvecm(i,:)),'_',int2str(ifig)])
% % end
% % end
% % close all
% end
end

View File

@ -55,7 +55,7 @@ else
load(filetoload,'lpmat','lpmat0','istable','iunstable','iindeterm','iwrong' ,'infox')
lpmat = [lpmat0 lpmat];
type = 'mc';
end
end
end
[Nsam, np] = size(lpmat);
npar = size(pnames,1);
@ -65,58 +65,91 @@ nbr_irf_restrictions = size(DynareOptions.endogenous_prior_restrictions.irf,1);
nbr_moment_restrictions = size(DynareOptions.endogenous_prior_restrictions.moment,1);
if init
mat_irf=cell(nbr_irf_restrictions,1);
for ij=1:nbr_irf_restrictions,
mat_irf{ij}=NaN(Nsam,length(DynareOptions.endogenous_prior_restrictions.irf{ij,3}));
end
mat_moment=cell(nbr_moment_restrictions,1);
for ij=1:nbr_moment_restrictions,
mat_moment{ij}=NaN(Nsam,length(DynareOptions.endogenous_prior_restrictions.moment{ij,3}));
end
irestrictions = [1:Nsam];
h = dyn_waitbar(0,'Please wait...');
for j=1:Nsam,
Model = set_all_parameters(lpmat(j,:)',EstimatedParameters,Model);
[Tt,Rr,SteadyState,info] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
if info(1)==0,
[info, info_irf, info_moment, data_irf, data_moment]=endogenous_prior_restrictions(Tt,Rr,Model,DynareOptions,DynareResults);
if ~isempty(info_irf)
for ij=1:nbr_irf_restrictions,
mat_irf{ij}(j,:)=data_irf{ij}(:,2)';
end
indx_irf(j,:)=info_irf(:,1);
end
if ~isempty(info_moment)
for ij=1:nbr_moment_restrictions,
mat_moment{ij}(j,:)=data_moment{ij}(:,2)';
end
indx_moment(j,:)=info_moment(:,1);
end
else
irestrictions(j)=0;
mat_irf=cell(nbr_irf_restrictions,1);
for ij=1:nbr_irf_restrictions,
mat_irf{ij}=NaN(Nsam,length(DynareOptions.endogenous_prior_restrictions.irf{ij,3}));
end
dyn_waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
end
dyn_waitbar_close(h);
irestrictions=irestrictions(find(irestrictions));
xmat=lpmat(irestrictions,:);
skipline()
endo_prior_restrictions=DynareOptions.endogenous_prior_restrictions;
save([OutputDirectoryName,filesep,fname_,'_',type,'_restrictions'],'xmat','mat_irf','mat_moment','irestrictions','indx_irf','indx_moment','endo_prior_restrictions');
else
load([OutputDirectoryName,filesep,fname_,'_',type,'_restrictions'],'xmat','mat_irf','mat_moment','irestrictions','indx_irf','indx_moment','endo_prior_restrictions');
end
if ~isempty(indx_irf),
% For single legend search which has maximum nbr of restrictions
mat_moment=cell(nbr_moment_restrictions,1);
for ij=1:nbr_moment_restrictions,
mat_moment{ij}=NaN(Nsam,length(DynareOptions.endogenous_prior_restrictions.moment{ij,3}));
end
irestrictions = [1:Nsam];
h = dyn_waitbar(0,'Please wait...');
for j=1:Nsam,
Model = set_all_parameters(lpmat(j,:)',EstimatedParameters,Model);
if nbr_moment_restrictions,
[Tt,Rr,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults);
else
[Tt,Rr,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
end
if info(1)==0,
[info, info_irf, info_moment, data_irf, data_moment]=endogenous_prior_restrictions(Tt,Rr,Model,DynareOptions,DynareResults);
if ~isempty(info_irf)
for ij=1:nbr_irf_restrictions,
mat_irf{ij}(j,:)=data_irf{ij}(:,2)';
end
indx_irf(j,:)=info_irf(:,1);
end
if ~isempty(info_moment)
for ij=1:nbr_moment_restrictions,
mat_moment{ij}(j,:)=data_moment{ij}(:,2)';
end
indx_moment(j,:)=info_moment(:,1);
end
else
irestrictions(j)=0;
end
dyn_waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
end
dyn_waitbar_close(h);
irestrictions=irestrictions(find(irestrictions));
xmat=lpmat(irestrictions,:);
skipline()
endo_prior_restrictions=DynareOptions.endogenous_prior_restrictions;
save([OutputDirectoryName,filesep,fname_,'_',type,'_restrictions'],'xmat','mat_irf','mat_moment','irestrictions','indx_irf','indx_moment','endo_prior_restrictions');
else
load([OutputDirectoryName,filesep,fname_,'_',type,'_restrictions'],'xmat','mat_irf','mat_moment','irestrictions','indx_irf','indx_moment','endo_prior_restrictions');
end
if ~isempty(indx_irf),
skipline()
disp('Deleting old IRF calibration plots ...')
a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_irf_calib*.eps']);
for j=1:length(a),
delete([OutputDirectoryName,filesep,a(j).name]);
end
a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_irf_calib*.fig']);
for j=1:length(a),
delete([OutputDirectoryName,filesep,a(j).name]);
end
a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_irf_calib*.pdf']);
for j=1:length(a),
delete([OutputDirectoryName,filesep,a(j).name]);
end
a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions.eps']);
for j=1:length(a),
delete([OutputDirectoryName,filesep,a(j).name]);
end
a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions.fig']);
for j=1:length(a),
delete([OutputDirectoryName,filesep,a(j).name]);
end
a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions.pdf']);
for j=1:length(a),
delete([OutputDirectoryName,filesep,a(j).name]);
end
disp('done !')
skipline()
% For single legend search which has maximum nbr of restrictions
all_irf_couples = cellstr([char(endo_prior_restrictions.irf(:,1)) char(endo_prior_restrictions.irf(:,2))]);
irf_couples = unique(all_irf_couples);
nbr_irf_couples = size(irf_couples,1);
plot_indx = NaN(nbr_irf_couples,1);
time_matrix=cell(nbr_irf_couples,1);
indx_irf_matrix=zeros(length(irestrictions),nbr_irf_couples);
irf_matrix=cell(nbr_irf_couples,1);
irf_mean=cell(nbr_irf_couples,1);
irf_median=cell(nbr_irf_couples,1);
@ -131,6 +164,7 @@ if ~isempty(indx_irf),
plot_indx(ij) = find(strcmp(irf_couples,all_irf_couples(ij,:)));
time_matrix{plot_indx(ij)} = [time_matrix{plot_indx(ij)} endo_prior_restrictions.irf{ij,3}];
end
iplot_indx = ones(size(plot_indx));
indx_irf = indx_irf(irestrictions,:);
h1=dyn_figure(DynareOptions,'name',[type ' evaluation of irf restrictions']);
@ -142,6 +176,7 @@ if ~isempty(indx_irf),
for ij=1:nbr_irf_restrictions,
mat_irf{ij}=mat_irf{ij}(irestrictions,:);
irf_matrix{plot_indx(ij)} = [irf_matrix{plot_indx(ij)} mat_irf{ij}];
indx_irf_matrix(:,plot_indx(ij)) = indx_irf_matrix(:,plot_indx(ij)) + indx_irf(:,ij);
for ik=1:size(mat_irf{ij},2),
[Mean,Median,Var,HPD,Distrib] = ...
posterior_moments(mat_irf{ij}(:,ik),0,DynareOptions.mh_conf_sig);
@ -156,7 +191,8 @@ if ~isempty(indx_irf),
if size(mat_irf{ij},2)>1,
leg = [leg,':' ,num2str(endo_prior_restrictions.irf{ij,3}(end))];
aleg = [aleg,'-' ,num2str(endo_prior_restrictions.irf{ij,3}(end))];
end
iplot_indx(ij)=0;
end
if length(time_matrix{plot_indx(ij)})==1,
figure(h1),
subplot(nrow,ncol, plot_indx(ij)),
@ -185,28 +221,28 @@ if ~isempty(indx_irf),
indx1 = find(indx_irf(:,ij)==0);
indx2 = find(indx_irf(:,ij)~=0);
atitle0=[endo_prior_restrictions.irf{ij,1},' vs ',endo_prior_restrictions.irf{ij,2}, '(', leg,')'];
fprintf(['%4.1f%% of the prior support matches IRF ',atitle0,' inside [%4.1f, %4.1f]\n'],length(indx1)/length(irestrictions)*100,endo_prior_restrictions.irf{ij,4})
fprintf(['%4.1f%% of the ',type,' support matches IRF ',atitle0,' inside [%4.1f, %4.1f]\n'],length(indx1)/length(irestrictions)*100,endo_prior_restrictions.irf{ij,4})
% aname=[type '_irf_calib_',int2str(ij)];
aname=[type '_irf_calib_',endo_prior_restrictions.irf{ij,1},'_vs_',endo_prior_restrictions.irf{ij,2},'_',aleg];
atitle=[type ' IRF Calib: Parameter(s) driving ',endo_prior_restrictions.irf{ij,1},' vs ',endo_prior_restrictions.irf{ij,2}, '(', leg,')'];
options_mcf.amcf_name = aname;
options_mcf.amcf_title = atitle;
options_mcf.beha_title = 'IRF prior restriction';
options_mcf.nobeha_title = 'NO IRF prior restriction';
options_mcf.beha_title = 'IRF restriction';
options_mcf.nobeha_title = 'NO IRF restriction';
options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2)
mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, DynareOptions);
end
% [proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
% indplot=find(proba<pvalue_ks);
% if ~isempty(indplot)
% stab_map_1(xmat, indx1, indx2, aname, 1, indplot, OutputDirectoryName,[],atitle);
% end
% [proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
% indplot=find(proba<pvalue_ks);
% if ~isempty(indplot)
% stab_map_1(xmat, indx1, indx2, aname, 1, indplot, OutputDirectoryName,[],atitle);
% end
end
figure(h1);
for ij=1:nbr_irf_couples,
if length(time_matrix{ij})>1,
figure(h1);
subplot(nrow,ncol, ij)
itmp = (find(plot_indx==ij));
plot(time_matrix{ij},[max(irf_matrix{ij})' min(irf_matrix{ij})'],'k--','linewidth',2)
@ -217,14 +253,14 @@ if ~isempty(indx_irf),
tmp=[];
for ir=1:length(itmp),
for it=1:length(endo_prior_restrictions.irf{itmp(ir),3})
temp_index = find(time_matrix{ij}==endo_prior_restrictions.irf{itmp(ir),3}(it));
tmp(temp_index,:) = endo_prior_restrictions.irf{itmp(ir),4};
temp_index = find(time_matrix{ij}==endo_prior_restrictions.irf{itmp(ir),3}(it));
tmp(temp_index,:) = endo_prior_restrictions.irf{itmp(ir),4};
end
end
% tmp = cell2mat(endo_prior_restrictions.irf(itmp,4));
% tmp = cell2mat(endo_prior_restrictions.irf(itmp,4));
tmp(isinf(tmp(:,1)),1)=a(3);
tmp(isinf(tmp(:,2)),2)=a(4);
hp = patch([time_matrix{ij} time_matrix{ij}(end:-1:1)],tmp(:),'b');
hp = patch([time_matrix{ij} time_matrix{ij}(end:-1:1)],[tmp(:,1); tmp(end:-1:1,2)],'b');
set(hp,'FaceAlpha',[0.5])
plot(a(1:2),[0 0],'r')
hold off,
@ -233,27 +269,79 @@ if ~isempty(indx_irf),
set(gca,'xtick',sort(time_matrix{ij}))
itmp = min(itmp);
title([endo_prior_restrictions.irf{itmp,1},' vs ',endo_prior_restrictions.irf{itmp,2}],'interpreter','none'),
if any(iplot_indx.*plot_indx==ij),
% MCF of the couples with logical AND
indx1 = find(indx_irf_matrix(:,ij)==0);
indx2 = find(indx_irf_matrix(:,ij)~=0);
leg = num2str(time_matrix{ij}(1));
leg = [leg '...' num2str(time_matrix{ij}(end))];
aleg = 'ALL';
atitle0=[endo_prior_restrictions.irf{itmp,1},' vs ',endo_prior_restrictions.irf{itmp,2}, '(', leg,')'];
fprintf(['%4.1f%% of the ',type,' support matches IRF restrictions ',atitle0,'\n'],length(indx1)/length(irestrictions)*100)
% aname=[type '_irf_calib_',int2str(ij)];
aname=[type '_irf_calib_',endo_prior_restrictions.irf{itmp,1},'_vs_',endo_prior_restrictions.irf{itmp,2},'_',aleg];
atitle=[type ' IRF Calib: Parameter(s) driving ',endo_prior_restrictions.irf{itmp,1},' vs ',endo_prior_restrictions.irf{itmp,2}, '(', leg,')'];
options_mcf.amcf_name = aname;
options_mcf.amcf_title = atitle;
options_mcf.beha_title = 'IRF restriction';
options_mcf.nobeha_title = 'NO IRF restriction';
options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2)
mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, DynareOptions);
end
end
end
end
dyn_saveas(h1,[OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions'],DynareOptions);
skipline()
end
if ~isempty(indx_moment)
skipline()
disp('Deleting old MOMENT calibration plots ...')
a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_moment_calib*.eps']);
for j=1:length(a),
delete([OutputDirectoryName,filesep,a(j).name]);
end
a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_moment_calib*.fig']);
for j=1:length(a),
delete([OutputDirectoryName,filesep,a(j).name]);
end
a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_moment_calib*.pdf']);
for j=1:length(a),
delete([OutputDirectoryName,filesep,a(j).name]);
end
a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions.eps']);
for j=1:length(a),
delete([OutputDirectoryName,filesep,a(j).name]);
end
a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions.fig']);
for j=1:length(a),
delete([OutputDirectoryName,filesep,a(j).name]);
end
a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions.pdf']);
for j=1:length(a),
delete([OutputDirectoryName,filesep,a(j).name]);
end
disp('done !')
skipline()
options_mcf.param_names = char(BayesInfo.name);
all_moment_couples = cellstr([char(endo_prior_restrictions.moment(:,1)) char(endo_prior_restrictions.moment(:,2))]);
moment_couples = unique(all_moment_couples);
nbr_moment_couples = size(moment_couples,1);
plot_indx = NaN(nbr_moment_couples,1);
time_matrix=cell(nbr_moment_couples,1);
indx_moment_matrix=zeros(length(irestrictions),nbr_moment_couples);
moment_matrix=cell(nbr_moment_couples,1);
moment_mean=cell(nbr_moment_couples,1);
moment_median=cell(nbr_moment_couples,1);
moment_var=cell(nbr_moment_couples,1);
moment_HPD=cell(nbr_moment_couples,1);
moment_distrib=cell(nbr_moment_couples,1);
% For single legend search which has maximum nbr of restrictions
% For single legend search which has maximum nbr of restrictions
maxijv=0;
for ij=1:nbr_moment_restrictions
if length(endo_prior_restrictions.moment{ij,3})>maxijv
@ -262,6 +350,7 @@ if ~isempty(indx_moment)
plot_indx(ij) = find(strcmp(moment_couples,all_moment_couples(ij,:)));
time_matrix{plot_indx(ij)} = [time_matrix{plot_indx(ij)} endo_prior_restrictions.moment{ij,3}];
end
iplot_indx = ones(size(plot_indx));
indx_moment = indx_moment(irestrictions,:);
h2=dyn_figure(DynareOptions,'name',[type ' evaluation of moment restrictions']);
@ -274,6 +363,7 @@ if ~isempty(indx_moment)
for ij=1:nbr_moment_restrictions,
mat_moment{ij}=mat_moment{ij}(irestrictions,:);
moment_matrix{plot_indx(ij)} = [moment_matrix{plot_indx(ij)} mat_moment{ij}];
indx_moment_matrix(:,plot_indx(ij)) = indx_moment_matrix(:,plot_indx(ij)) + indx_moment(:,ij);
for ik=1:size(mat_moment{ij},2),
[Mean,Median,Var,HPD,Distrib] = ...
posterior_moments(mat_moment{ij}(:,ik),0,DynareOptions.mh_conf_sig);
@ -288,6 +378,7 @@ if ~isempty(indx_moment)
if size(mat_moment{ij},2)>1,
leg = [leg,':' ,num2str(endo_prior_restrictions.moment{ij,3}(end))];
aleg = [aleg,'_' ,num2str(endo_prior_restrictions.moment{ij,3}(end))];
iplot_indx(ij)=0;
end
if length(time_matrix{plot_indx(ij)})==1,
figure(h2),
@ -303,38 +394,38 @@ if ~isempty(indx_moment)
set(hp,'FaceAlpha', 0.5)
hold off,
title([endo_prior_restrictions.moment{ij,1},' vs ',endo_prior_restrictions.moment{ij,2},'(',leg,')'],'interpreter','none'),
% if ij==maxij
% leg1 = num2str(endo_prior_restrictions.moment{ij,3}(:));
% [legend_h,object_h,plot_h,text_strings]=legend(leg1);
% Position=get(legend_h,'Position');Position(1:2)=[-0.055 0.95-Position(4)];
% set(legend_h,'Position',Position);
% end
% if ij==maxij
% leg1 = num2str(endo_prior_restrictions.moment{ij,3}(:));
% [legend_h,object_h,plot_h,text_strings]=legend(leg1);
% Position=get(legend_h,'Position');Position(1:2)=[-0.055 0.95-Position(4)];
% set(legend_h,'Position',Position);
% end
end
indx1 = find(indx_moment(:,ij)==0);
indx2 = find(indx_moment(:,ij)~=0);
atitle0=[endo_prior_restrictions.moment{ij,1},' vs ',endo_prior_restrictions.moment{ij,2}, '(', leg,')'];
fprintf(['%4.1f%% of the prior support matches MOMENT ',atitle0,' inside [%4.1f, %4.1f]\n'],length(indx1)/length(irestrictions)*100,endo_prior_restrictions.moment{ij,4})
fprintf(['%4.1f%% of the ',type,' support matches MOMENT ',atitle0,' inside [%4.1f, %4.1f]\n'],length(indx1)/length(irestrictions)*100,endo_prior_restrictions.moment{ij,4})
% aname=[type '_moment_calib_',int2str(ij)];
aname=[type '_moment_calib_',endo_prior_restrictions.moment{ij,1},'_vs_',endo_prior_restrictions.moment{ij,2},'_',aleg];
atitle=[type ' MOMENT Calib: Parameter(s) driving ',endo_prior_restrictions.moment{ij,1},' vs ',endo_prior_restrictions.moment{ij,2}, '(', leg,')'];
options_mcf.amcf_name = aname;
options_mcf.amcf_title = atitle;
options_mcf.beha_title = 'moment prior restriction';
options_mcf.nobeha_title = 'NO moment prior restriction';
options_mcf.beha_title = 'moment restriction';
options_mcf.nobeha_title = 'NO moment restriction';
options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2)
mcf_analysis(xmat, indx1, indx2, options_mcf, DynareOptions);
end
% [proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
% indplot=find(proba<pvalue_ks);
% if ~isempty(indplot)
% stab_map_1(xmat, indx1, indx2, aname, 1, indplot, OutputDirectoryName,[],atitle);
% end
% [proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
% indplot=find(proba<pvalue_ks);
% if ~isempty(indplot)
% stab_map_1(xmat, indx1, indx2, aname, 1, indplot, OutputDirectoryName,[],atitle);
% end
end
figure(h2);
for ij=1:nbr_moment_couples,
if length(time_matrix{ij})>1,
figure(h2);
subplot(nrow,ncol, ij)
itmp = (find(plot_indx==ij));
plot(time_matrix{ij},[max(moment_matrix{ij})' min(moment_matrix{ij})'],'k--','linewidth',2)
@ -345,14 +436,14 @@ if ~isempty(indx_moment)
tmp=[];
for ir=1:length(itmp),
for it=1:length(endo_prior_restrictions.moment{itmp(ir),3})
temp_index = find(time_matrix{ij}==endo_prior_restrictions.moment{itmp(ir),3}(it));
tmp(temp_index,:) = endo_prior_restrictions.moment{itmp(ir),4};
temp_index = find(time_matrix{ij}==endo_prior_restrictions.moment{itmp(ir),3}(it));
tmp(temp_index,:) = endo_prior_restrictions.moment{itmp(ir),4};
end
end
% tmp = cell2mat(endo_prior_restrictions.moment(itmp,4));
% tmp = cell2mat(endo_prior_restrictions.moment(itmp,4));
tmp(isinf(tmp(:,1)),1)=a(3);
tmp(isinf(tmp(:,2)),2)=a(4);
hp = patch([time_matrix{ij} time_matrix{ij}(end:-1:1)],tmp(:),'b');
hp = patch([time_matrix{ij} time_matrix{ij}(end:-1:1)],[tmp(:,1); tmp(end:-1:1,2)],'b');
set(hp,'FaceAlpha',[0.5])
plot(a(1:2),[0 0],'r')
hold off,
@ -361,6 +452,27 @@ if ~isempty(indx_moment)
set(gca,'xtick',sort(time_matrix{ij}))
itmp = min(itmp);
title([endo_prior_restrictions.moment{itmp,1},' vs ',endo_prior_restrictions.moment{itmp,2}],'interpreter','none'),
if any(iplot_indx.*plot_indx==ij),
% MCF of the couples with logical AND
indx1 = find(indx_moment_matrix(:,ij)==0);
indx2 = find(indx_moment_matrix(:,ij)~=0);
leg = num2str(time_matrix{ij}(1));
leg = [leg '...' num2str(time_matrix{ij}(end))];
aleg = 'ALL';
atitle0=[endo_prior_restrictions.moment{itmp,1},' vs ',endo_prior_restrictions.moment{itmp,2}, '(', leg,')'];
fprintf(['%4.1f%% of the ',type,' support matches MOMENT restrictions ',atitle0,'\n'],length(indx1)/length(irestrictions)*100)
% aname=[type '_moment_calib_',int2str(ij)];
aname=[type '_moment_calib_',endo_prior_restrictions.moment{itmp,1},'_vs_',endo_prior_restrictions.moment{itmp,2},'_',aleg];
atitle=[type ' MOMENT Calib: Parameter(s) driving ',endo_prior_restrictions.moment{itmp,1},' vs ',endo_prior_restrictions.moment{itmp,2}, '(', leg,')'];
options_mcf.amcf_name = aname;
options_mcf.amcf_title = atitle;
options_mcf.beha_title = 'moment restriction';
options_mcf.nobeha_title = 'NO moment restriction';
options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2)
mcf_analysis(xmat, indx1, indx2, options_mcf, DynareOptions);
end
end
end
end
dyn_saveas(h2,[OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions'],DynareOptions);

View File

@ -267,7 +267,11 @@ if fload==0,
M_ = set_all_parameters([lpmat0(j,:) lpmat(j,:)]',estim_params_,M_);
%try stoch_simul([]);
try
[Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
if ~ isempty(options_.endogenous_prior_restrictions.moment)
[Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
else
[Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
end
infox(j,1)=info(1);
if infox(j,1)==0 && ~exist('T'),
dr_=oo_.dr;

View File

@ -273,6 +273,8 @@ for b=fpar:B
if irun(5) > MAX_nruns || b == B
stock = stock_param(1:irun(5)-1,:);
stock_logpo = stock_logpo(1:irun(5)-1);
stock_ys = stock_ys(1:irun(5)-1,:);
ifil(5) = ifil(5) + 1;
save([DirectoryName '/' M_.fname '_param' int2str(ifil(5)) '.mat'],'stock','stock_logpo','stock_ys');
if RemoteFlag==1,

View File

@ -7,6 +7,7 @@ MODFILES = \
estimation/fs2000_initialize_from_calib.mod \
estimation/fs2000_calibrated_covariance.mod \
gsa/ls2003.mod \
gsa/ls2003a.mod \
ramst.mod \
ramst_a.mod \
ramst_static_tag.mod \

View File

@ -54,11 +54,18 @@ rho_q ,beta_pdf,0.4,0.2;
rho_A ,beta_pdf,0.5,0.2;
rho_ys ,beta_pdf,0.8,0.1;
rho_pies,beta_pdf,0.7,0.15;
stderr e_R,inv_gamma_pdf,1.2533,0.6551;
stderr e_q,inv_gamma_pdf,2.5066,1.3103;
stderr e_A,inv_gamma_pdf,1.2533,0.6551;
stderr e_ys,inv_gamma_pdf,1.2533,0.6551;
stderr e_pies,inv_gamma_pdf,1.88,0.9827;
/*
stderr e_R,inv_gamma_pdf,(1.2533),(0.6551);
stderr e_q,inv_gamma_pdf,(2.5066),(1.3103);
stderr e_A,inv_gamma_pdf,(1.2533),(0.6551);
stderr e_ys,inv_gamma_pdf,(1.2533),(0.6551);
stderr e_pies,inv_gamma_pdf,(1.88),(0.9827);
*/
stderr e_R,inv_gamma_pdf,(1.2533/3),(0.6551/10);
stderr e_q,inv_gamma_pdf,(2.5066/3),(1.3103/10);
stderr e_A,inv_gamma_pdf,(1.2533/3),(0.6551/10);
stderr e_ys,inv_gamma_pdf,(1.2533/3),(0.6551/10);
stderr e_pies,inv_gamma_pdf,(1.88/3),(0.9827/10);
end;
// endogenous prior restrictions
@ -83,12 +90,16 @@ end;
moment_calibration;
//y_obs,y_obs, [0.8 1.1]; //[unconditional variance]
y_obs,y_obs(-(1:4)), +; //[first year acf]
@#for ilag in 0:1
y_obs,R_obs(-@{ilag}), -; //[ccf]
//y_obs,pie_obs(-4:4), -; //[ccf]
@#for ilag in -2:2
y_obs,R_obs(@{ilag}), -; //[ccf]
@#endfor
@#for ilag in -4:4
y_obs,pie_obs(@{ilag}), -; //[ccf]
@#endfor
end;
dynare_sensitivity(prior_range=0);
dynare_sensitivity(prior_range=0, nodisplay, graph_format=(fig,eps));
/*
estimation(datafile='data_ca1.m',first_obs=8,nobs=79,mh_nblocks=2, mode_file = ls2003a_mode,
prefilter=1,mh_jscale=0.5,mh_replic=5000, mode_compute=0, mh_drop=0.6, bayesian_irf);