diff --git a/matlab/dynare_sensitivity.m b/matlab/dynare_sensitivity.m
index d85ca1f22bff00ac3d04c77f88535c43101d1063..86c2a8caa0f47b3cce845eaafdb148cffb5be11b 100644
--- a/matlab/dynare_sensitivity.m
+++ b/matlab/dynare_sensitivity.m
@@ -145,13 +145,13 @@ options_gsa = set_default_option(options_gsa,'load_stab',0);
options_gsa = set_default_option(options_gsa,'alpha2_stab',0);
options_gsa = set_default_option(options_gsa,'ksstat',0.1);
options_gsa = set_default_option(options_gsa,'pvalue_ks',0.001);
-options_gsa = set_default_option(options_gsa,'pvalue_corr',0.001);
+options_gsa = set_default_option(options_gsa,'pvalue_corr',1.e-5);
%options_gsa = set_default_option(options_gsa,'load_mh',0);
% REDFORM mapping
options_gsa = set_default_option(options_gsa,'logtrans_redform',0);
options_gsa = set_default_option(options_gsa,'threshold_redform',[]);
options_gsa = set_default_option(options_gsa,'ksstat_redform',0.001);
-options_gsa = set_default_option(options_gsa,'alpha2_redform',0);
+options_gsa = set_default_option(options_gsa,'alpha2_redform',1.e-5);
options_gsa = set_default_option(options_gsa,'namendo',[]);
options_gsa = set_default_option(options_gsa,'namlagendo',[]);
options_gsa = set_default_option(options_gsa,'namexo',[]);
@@ -161,7 +161,7 @@ options_gsa = set_default_option(options_gsa,'var_rmse',char(options_.varobs));
options_gsa = set_default_option(options_gsa,'pfilt_rmse',0.1);
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',0);
+options_gsa = set_default_option(options_gsa,'alpha2_rmse',1.e-5);
if options_gsa.redform && options_gsa.neighborhood_width==0 && isempty(options_gsa.threshold_redform),
options_gsa.pprior=1;
diff --git a/matlab/gsa/map_calibration.m b/matlab/gsa/map_calibration.m
index 5f9df099daefef5f2058665800b54f9b2e47df46..fdd0ab142ab48572d5123b82433d6c2e0d33a653 100644
--- a/matlab/gsa/map_calibration.m
+++ b/matlab/gsa/map_calibration.m
@@ -22,7 +22,7 @@ pnames = Model.param_names(EstimatedParameters.param_vals(:,1),:);
pvalue_ks = DynareOptions.opt_gsa.pvalue_ks;
indx_irf = [];
indx_moment = [];
-
+DynareOptions.nodisplay = 1;
skipline()
disp('Sensitivity analysis for calibration criteria')
@@ -63,6 +63,7 @@ for ij=1:nbr_moment_restrictions,
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');
@@ -83,11 +84,15 @@ for j=1:Nsam,
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()
-save([OutputDirectoryName,filesep,fname_,'_',type,'_restrictions'],'xmat','mat_irf','mat_moment','irestrictions','indx_irf','indx_moment');
+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');
if ~isempty(indx_irf),
indx_irf = indx_irf(irestrictions,:);
@@ -97,6 +102,13 @@ if ~isempty(indx_irf),
if nrow*(nrow-1)>nbr_irf_restrictions,
ncol=nrow-1;
end
+ % For single legend search which has maximum nbr of restrictions
+ maxijv=0;
+ for ij=1:nbr_irf_restrictions
+ if length(DynareOptions.endogenous_prior_restrictions.irf{ij,3})>maxijv
+ maxij=ij;maxijv=length(DynareOptions.endogenous_prior_restrictions.irf{ij,3});
+ end
+ end
for ij=1:nbr_irf_restrictions,
figure(h1),
mat_irf{ij}=mat_irf{ij}(irestrictions,:);
@@ -115,14 +127,24 @@ if ~isempty(indx_irf),
leg = num2str(DynareOptions.endogenous_prior_restrictions.irf{ij,3}(1));
if size(mat_irf{ij},2)>1,
leg = [leg,':' ,num2str(DynareOptions.endogenous_prior_restrictions.irf{ij,3}(end))];
- end
- title([DynareOptions.endogenous_prior_restrictions.irf{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.irf{ij,2}, '(', leg,')'],'interpreter','none'),
-
+ end
+ title([DynareOptions.endogenous_prior_restrictions.irf{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.irf{ij,2}, '(', leg,')'],'interpreter','none'),
+ %set(legend_h,'Xlim',[0 1]);
+ if ij==maxij
+ leg1 = num2str(DynareOptions.endogenous_prior_restrictions.irf{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
+ % hc = get(h,'Children');
+ %for i=2:2:length(hc)
+ %end
indx1 = find(indx_irf(:,ij)==0);
indx2 = find(indx_irf(:,ij)~=0);
atitle=[DynareOptions.endogenous_prior_restrictions.irf{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.irf{ij,2}, '(', leg,')'];
fprintf(['%4.1f%% of the prior support matches IRF ',atitle,' inside [%4.1f, %4.1f]\n'],length(indx1)/length(irestrictions)*100,DynareOptions.endogenous_prior_restrictions.irf{ij,4})
- aname=[type '_irf_calib_',int2str(ij)];
+ % aname=[type '_irf_calib_',int2str(ij)];
+ aname=[type '_irf_calib_',endo_prior_restrictions.irf{ij,1},'_VS_',endo_prior_restrictions.irf{ij,2}];
atitle=[type ' IRF Calib: Parameter(s) driving ',DynareOptions.endogenous_prior_restrictions.irf{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.irf{ij,2}, '(', leg,')'];
[proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
indplot=find(proba<pvalue_ks);
@@ -142,6 +164,13 @@ if ~isempty(indx_moment)
if nrow*(nrow-1)>nbr_moment_restrictions,
ncol=nrow-1;
end
+ % For single legend search which has maximum nbr of restrictions
+ maxijv=0;
+ for ij=1:nbr_moment_restrictions
+ if length(DynareOptions.endogenous_prior_restrictions.moment{ij,3})>maxijv
+ maxij=ij;maxijv=length(DynareOptions.endogenous_prior_restrictions.moment{ij,3});
+ end
+ end
for ij=1:nbr_moment_restrictions,
figure(h2),
mat_moment{ij}=mat_moment{ij}(irestrictions,:);
@@ -162,12 +191,18 @@ if ~isempty(indx_moment)
leg = [leg,':' ,num2str(DynareOptions.endogenous_prior_restrictions.moment{ij,3}(end))];
end
title([DynareOptions.endogenous_prior_restrictions.moment{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.moment{ij,2},'(',leg,')'],'interpreter','none'),
-
+ if ij==maxij
+ leg1 = num2str(DynareOptions.endogenous_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
indx1 = find(indx_moment(:,ij)==0);
indx2 = find(indx_moment(:,ij)~=0);
atitle=[DynareOptions.endogenous_prior_restrictions.moment{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.moment{ij,2}, '(', leg,')'];
fprintf(['%4.1f%% of the prior support matches MOMENT ',atitle,' inside [%4.1f, %4.1f]\n'],length(indx1)/length(irestrictions)*100,DynareOptions.endogenous_prior_restrictions.moment{ij,4})
- aname=[type '_moment_calib_',int2str(ij)];
+ % aname=[type '_moment_calib_',int2str(ij)];
+ aname=[type '_moment_calib_',DynareOptions.endogenous_prior_restrictions.moment(ij,1),'_VS_',DynareOptions.endogenous_prior_restrictions.moment(ij,2)];
atitle=[type ' MOMENT Calib: Parameter(s) driving ',DynareOptions.endogenous_prior_restrictions.moment{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.moment{ij,2}, '(', leg,')'];
[proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
indplot=find(proba<pvalue_ks);
diff --git a/matlab/gsa/scatter_mcf.m b/matlab/gsa/scatter_mcf.m
new file mode 100644
index 0000000000000000000000000000000000000000..b753a374e684998f64bb86d44dcbfbc3b59cbc18
--- /dev/null
+++ b/matlab/gsa/scatter_mcf.m
@@ -0,0 +1,163 @@
+function scatter_mcf(X,Y,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, DynareOptions, beha_name, non_beha_name)
+% Frontend to the Sensitivity Analysis Toolbox for DYNARE
+%
+% Reference:
+% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
+
+% Copyright (C) 2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+
+% PURPOSE: Pairwise scatter plots of the columns of x and y after
+% Monte Carlo filtering
+%---------------------------------------------------
+% USAGE: scatter_mcf(x,y,vnames,pltsym,diagon)
+% or scatter_mcf(x,y) which relies on defaults
+% where:
+% x = an nxk matrix with columns containing behavioural sample
+% y = an mxk matrix with columns containing non-behavioural sample
+% vnames = a vector of variable names
+% (default = numeric labels 1,2,3 etc.)
+% pltsym = a plt symbol
+% (default = '.' for npts > 100, 'o' for npts < 100
+
+
+Z=[X;Y];
+[n,p] = size(X);
+% X = X - ones(n,1)*min(Z);
+% X = X ./ (ones(n,1)*max(Z));
+[n,p] = size(Y);
+% Y = Y - ones(n,1)*min(Z);
+% Y = Y ./ (ones(n,1)*max(Z));
+[n,p] = size(Z);
+clear Z;
+
+nflag = 0;
+if nargin >=3
+ nflag = 1;
+end;
+
+if nargin<4 || isempty(plotsymbol)
+ if n*p<100, plotsymbol = 'o';
+ else plotsymbol = '.';
+ end
+end
+
+if nargin<5
+ fnam='';
+end
+if nargin<6,
+ dirname='';
+ nograph=1;
+else
+ nograph=0;
+end
+if nargin<7,
+ figtitle=fnam;
+end
+if nargin<8,
+ xparam1=[];
+end
+if nargin<10,
+ beha_name = 'BEHAVIOUR';
+ non_beha_name = 'NON-BEHAVIOUR';
+end
+if nargin==10,
+ non_beha_name = ['NON-' beha_name];
+end
+
+fig_nam_=[fnam];
+if ~nograph,
+ hh=dyn_figure(DynareOptions,'name',figtitle);
+end
+
+bf = 0.1;
+ffs = 0.05/(p-1);
+ffl = (1-2*bf-0.05)/p;
+fL = linspace(bf,1-bf+ffs,p+1);
+for i = 1:p
+ for j = 1:p
+ h = axes('position',[fL(i),fL(p+1-j),ffl,ffl]);
+ if i==j
+ h1=cumplot(X(:,j));
+ % set(h1,'color',[0 0 1], 'linestyle','--','LineWidth',1.5)
+ set(h1,'color',[0 0 1],'LineWidth',1.5)
+ hold on,
+ h2=cumplot(Y(:,j));
+ set(h2,'color',[1 0 0],'LineWidth',1.5)
+ if ~isempty(xparam1)
+ hold on, plot(xparam1([j j]),[0 1],'k--')
+ end
+ if j<p
+ set(gca,'XTickLabel',[],'XTick',[]);
+ else
+ grid off
+ end
+ set(gca,'YTickLabel',[],'YTick',[]);
+ else
+ if j>i
+ plot(X(:,i),X(:,j),[plotsymbol,'b'])
+ hold on,
+ plot(Y(:,i),Y(:,j),[plotsymbol,'r'])
+ else
+ plot(Y(:,i),Y(:,j),[plotsymbol,'r'])
+ hold on,
+ plot(X(:,i),X(:,j),[plotsymbol,'b'])
+ end
+ if ~isempty(xparam1)
+ hold on, plot(xparam1(i),xparam1(j),'k*')
+ end
+ hold off;
+ % axis([-0.1 1.1 -0.1 1.1])
+ if i<p,
+ set(gca,'YTickLabel',[],'YTick',[]);
+ else
+ set(gca,'yaxislocation','right');
+ end
+ if j<p
+ set(gca,'XTickLabel',[],'XTick',[]);
+ end
+ end
+ if nflag == 1
+ set(gca,'fontsize',9);
+ end;
+ if i==1
+ if nflag == 1
+ ylabel(vnames(j,:),'Rotation',45,'interpreter','none', ...
+ 'HorizontalAlignment','right','VerticalAlignment','middle');
+ else
+ ylabel([num2str(j),' '],'Rotation',90)
+ end;
+ end
+ if j==1
+ if nflag == 1
+ title(vnames(i,:),'interpreter','none','Rotation',45, ...
+ 'HorizontalAlignment','left','VerticalAlignment','bottom')
+ else
+ title(num2str(i))
+ end;
+ end
+ drawnow
+ end
+end
+if ~isoctave
+ annotation('textbox', [0.1,0,0.35,0.05],'String', beha_name,'Color','Blue','horizontalalignment','center');
+ annotation('textbox', [0.55,0,0.35,0.05],'String', non_beha_name,'Color','Red','horizontalalignment','center');
+end
+
+if ~nograph,
+ dyn_saveas(hh,[dirname,filesep,fig_nam_],DynareOptions);
+end
\ No newline at end of file
diff --git a/matlab/gsa/stab_map_.m b/matlab/gsa/stab_map_.m
index 92016f031d9a0a9de7b188605cd12f1282b7fb79..1b27045744b0b4e141a653b771e92295dd820d17 100644
--- a/matlab/gsa/stab_map_.m
+++ b/matlab/gsa/stab_map_.m
@@ -254,14 +254,13 @@ if fload==0,
iunstable=[1:Nsam];
iindeterm=zeros(1,Nsam);
iwrong=zeros(1,Nsam);
+ inorestriction=zeros(1,Nsam);
+ irestriction=zeros(1,Nsam);
for j=1:Nsam,
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 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;
@@ -303,7 +302,6 @@ if fload==0,
dr_ = oo_.dr;
if isfield(dr_,'ghx'),
egg(:,j) = sort(dr_.eigval);
- iunstable(j)=0;
if prepSA
jstab=jstab+1;
T(:,:,jstab) = [dr_.ghx dr_.ghu];
@@ -317,6 +315,15 @@ if fload==0,
nboth = dr_.nboth;
nfwrd = dr_.nfwrd;
end
+ info=endogenous_prior_restrictions(Tt,Rr,M_,options_,oo_);
+ infox(j,1)=info(1);
+ if info(1),
+ iwrong(j)=j;
+ inorestriction(j)=j;
+ else
+ iunstable(j)=0;
+ irestriction(j)=j;
+ end
else
istable(j)=0;
if isfield(dr_,'eigval')
@@ -348,10 +355,13 @@ if fload==0,
else
T=[];
end
- istable=istable(find(istable)); % stable params
- iunstable=iunstable(find(iunstable)); % unstable params
+ istable=istable(find(istable)); % stable params ignoring restrictions
+ irestriction=irestriction(find(irestriction)); % stable params & restrictions OK
+ inorestriction=inorestriction(find(inorestriction)); % stable params violating restrictions
+ iunstable=iunstable(find(iunstable)); % violation of BK & restrictions & solution could not be found (whatever goes wrong)
iindeterm=iindeterm(find(iindeterm)); % indeterminacy
iwrong=iwrong(find(iwrong)); % dynare could not find solution
+ ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong]))); % explosive roots
% % map stable samples
% istable=[1:Nsam];
@@ -395,22 +405,22 @@ if fload==0,
if pprior,
if ~prepSA
save([OutputDirectoryName filesep fname_ '_prior.mat'], ...
- 'bkpprior','lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
+ 'bkpprior','lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun', ...
'egg','yys','nspred','nboth','nfwrd','infox')
else
save([OutputDirectoryName filesep fname_ '_prior.mat'], ...
- 'bkpprior','lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
+ 'bkpprior','lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun', ...
'egg','yys','T','nspred','nboth','nfwrd','infox')
end
else
if ~prepSA
save([OutputDirectoryName filesep fname_ '_mc.mat'], ...
- 'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
+ 'lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun', ...
'egg','yys','nspred','nboth','nfwrd','infox')
else
save([OutputDirectoryName filesep fname_ '_mc.mat'], ...
- 'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
+ 'lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun', ...
'egg','yys','T','nspred','nboth','nfwrd','infox')
end
end
@@ -420,7 +430,7 @@ else
else
filetoload=[OutputDirectoryName filesep fname_ '_mc.mat'];
end
- load(filetoload,'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','nspred','nboth','nfwrd','infox')
+ load(filetoload,'lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun','egg','yys','nspred','nboth','nfwrd','infox')
Nsam = size(lpmat,1);
if pprior==0,
eval(['load ' options_.mode_file '.mat;']);
@@ -472,11 +482,13 @@ if pprior
aindetname=[aname, '_indet']; aindettitle='Prior StabMap: Parameter driving indeterminacy';
aunstablename=[aname, '_unst']; aunstabletitle='Prior StabMap: Parameter driving explosiveness of solution';
awronguniname=[aname, '_wrong']; awrongunititle='Prior StabMap: Parameter driving inability to find solution';
+ acalibuniname=[aname, '_calib']; acalibunititle='Prior StabMap: Parameter driving IRF/moment restrictions';
% bivariate
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';
+ acalibname='prior_calib'; acalibtitle='Prior StabMap: IRF/moment restrictions';
asname='prior_stable'; astitle='Prior StabMap: unique Stable Saddle-Path';
else
% univariate
@@ -484,11 +496,13 @@ else
aindetname=[aname, '_indet']; aindettitle='MC (around posterior mode) StabMap: Parameter driving indeterminacy';
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';
+ acalibuniname=[aname, '_calib']; acalibunititle='MC (around posterior mode) StabMap: Parameter driving IRF/moment restrictions';
% bivariate
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';
+ acalibname='mc_calib'; acalibtitle='MC (around posterior mode) StabMap: IRF/moment restrictions';
asname='mc_stable'; astitle='MC (around posterior mode) StabMap: Unique Stable Saddle-Path';
end
delete([OutputDirectoryName,filesep,fname_,'_',aname,'_*.*']);
@@ -498,7 +512,7 @@ delete([OutputDirectoryName,filesep,fname_,'_',auname,'_corr_*.*']);
delete([OutputDirectoryName,filesep,fname_,'_',aunstname,'_corr_*.*']);
delete([OutputDirectoryName,filesep,fname_,'_',aindname,'_corr_*.*']);
-if length(iunstable)>0 ,
+if length(iunstable)>0 || length(iwrong)>0,
fprintf(['%4.1f%% of the prior support gives unique saddle-path solution.\n'],length(istable)/Nsam*100)
fprintf(['%4.1f%% of the prior support gives explosive dynamics.\n'],(length(iunstable)-length(iwrong)-length(iindeterm) )/Nsam*100)
if ~isempty(iindeterm),
@@ -544,99 +558,182 @@ if length(iunstable)>0 ,
end
skipline()
- if length(iunstable)<Nsam
+ if length(iunstable)<Nsam || length(istable)>1
% Blanchard Kahn
- [proba, dproba] = stab_map_1(lpmat, istable, iunstable, aname,0);
+ itot = [1:Nsam];
+ itmp = itot(find(~ismember(itot,istable)));
+ [proba, dproba] = stab_map_1(lpmat, istable, itmp, aname,0);
% indstab=find(dproba>ksstat);
indstab=find(proba<pvalue_ks);
- disp('Smirnov statistics in driving acceptable behaviour')
+ [tmp,jtmp] = sort(proba(indstab),2,'ascend');
+ indstab = indstab(jtmp);
+ disp('Smirnov statistics in driving unique solution')
for j=1:length(indstab),
disp([M_.param_names(estim_params_.param_vals(indstab(j),1),:),' d-stat = ', num2str(dproba(indstab(j)),'%1.3f'),' p-value = ', num2str(proba(indstab(j)),'%1.3f')])
end
skipline()
- if ~isempty(indstab)
- stab_map_1(lpmat, istable, iunstable, aname, 1, indstab, OutputDirectoryName,[],atitle);
+ indcorr1 = stab_map_2(lpmat(istable,:),alpha2, pvalue_corr, asname);
+ indcorr2 = stab_map_2(lpmat(itmp,:),alpha2, pvalue_corr, auname);
+ indcorr = union(indcorr1(:), indcorr2(:));
+ indcorr = indcorr(~ismember(indcorr(:),indstab));
+ indstab = [indstab(:); indcorr(:)];
+ skipline()
+ if ~isempty(indstab) && ~options_.nograph
+% stab_map_1(lpmat, istable, iunstable, aname, 1, indstab, OutputDirectoryName,[],atitle);
+ scatter_mcf(lpmat(istable,indstab),lpmat(itmp,indstab),cellstr(bayestopt_.name(indstab+nshock)), ...
+ '.', [fname_,'_',asname], OutputDirectoryName,atitle,[], options_, ...
+ 'unique Stable Saddle-Path', 'NO unique Stable Saddle-Path')
end
- ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong])));
+% ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong])));
if ~isempty(iindeterm),
- [proba, dproba] = stab_map_1(lpmat, [1:Nsam], iindeterm, aindetname ,0);
+ itmp = itot(find(~ismember(itot,iindeterm)));
+ [proba, dproba] = stab_map_1(lpmat, itmp, iindeterm, aindetname ,0);
% indindet=find(dproba>ksstat);
indindet=find(proba<pvalue_ks);
+ [tmp,jtmp] = sort(proba(indindet),2,'ascend');
+ indindet = indindet(jtmp);
disp('Smirnov statistics in driving indeterminacy')
for j=1:length(indindet),
disp([M_.param_names(estim_params_.param_vals(indindet(j),1),:),' d-stat = ', num2str(dproba(indindet(j)),'%1.3f'),' p-value = ', num2str(proba(indindet(j)),'%1.3f')])
end
skipline()
+ indcorr1 = stab_map_2(lpmat(iindeterm,:),alpha2, pvalue_corr, aindname);
+ indcorr2 = stab_map_2(lpmat(itmp,:),alpha2, pvalue_corr, ['NOT_',aindname]);
+ indcorr = union(indcorr1(:), indcorr2(:));
+ indcorr = indcorr(~ismember(indcorr(:),indindet));
+ indindet = [indindet(:); indcorr(:)];
+ skipline()
if ~isempty(indindet)
- stab_map_1(lpmat, [1:Nsam], iindeterm, aindetname, 1, indindet, OutputDirectoryName,[],aindettitle);
+% stab_map_1(lpmat, itmp, iindeterm, aindetname, 1, indindet, OutputDirectoryName,[],aindettitle);
+ scatter_mcf(lpmat(itmp,indindet),lpmat(iindeterm,indindet),cellstr(bayestopt_.name(indindet+nshock)), ...
+ '.', [fname_,'_',aindname], OutputDirectoryName,aindtitle,[], options_, ...
+ 'NO indeterminacy', 'indeterminacy')
end
end
if ~isempty(ixun),
- [proba, dproba] = stab_map_1(lpmat, [1:Nsam], ixun, aunstablename,0);
+ itmp = itot(find(~ismember(itot,ixun)));
+ [proba, dproba] = stab_map_1(lpmat, itmp, ixun, aunstablename,0);
% indunst=find(dproba>ksstat);
indunst=find(proba<pvalue_ks);
+ [tmp,jtmp] = sort(proba(indunst),2,'ascend');
+ indunst = indunst(jtmp);
disp('Smirnov statistics in driving instability')
for j=1:length(indunst),
disp([M_.param_names(estim_params_.param_vals(indunst(j),1),:),' d-stat = ', num2str(dproba(indunst(j)),'%1.3f'),' p-value = ', num2str(proba(indunst(j)),'%1.3f')])
end
skipline()
+ indcorr1 = stab_map_2(lpmat(ixun,:),alpha2, pvalue_corr, aunstname);
+ indcorr2 = stab_map_2(lpmat(itmp,:),alpha2, pvalue_corr, ['NOT_',aunstname]);
+ indcorr = union(indcorr1(:), indcorr2(:));
+ indcorr = indcorr(~ismember(indcorr(:),indunst));
+ indunst = [indunst(:); indcorr(:)];
+ skipline()
if ~isempty(indunst)
- stab_map_1(lpmat, [1:Nsam], ixun, aunstablename, 1, indunst, OutputDirectoryName,[],aunstabletitle);
+% stab_map_1(lpmat, itmp, ixun, aunstablename, 1, indunst, OutputDirectoryName,[],aunstabletitle);
+ scatter_mcf(lpmat(itmp,indunst),lpmat(ixun,indunst),cellstr(bayestopt_.name(indunst+nshock)), ...
+ '.', [fname_,'_',aunstname], OutputDirectoryName,aunsttitle,[], options_, ...
+ 'NO explosive solution', 'explosive solution')
end
end
+ inorestriction = istable(find(~ismember(istable,irestriction))); % what went wrong beyong prior restrictions
+ iwrong = iwrong(find(~ismember(iwrong,inorestriction))); % what went wrong beyong prior restrictions
if ~isempty(iwrong),
- [proba, dproba] = stab_map_1(lpmat, [1:Nsam], iwrong, awronguniname,0);
+ itmp = itot(find(~ismember(itot,iwrong)));
+ [proba, dproba] = stab_map_1(lpmat, itmp, iwrong, awronguniname,0);
% indwrong=find(dproba>ksstat);
indwrong=find(proba<pvalue_ks);
+ [tmp,jtmp] = sort(proba(indwrong),2,'ascend');
+ indwrong = indwrong(jtmp);
disp('Smirnov statistics in driving no solution')
for j=1:length(indwrong),
disp([M_.param_names(estim_params_.param_vals(indwrong(j),1),:),' d-stat = ', num2str(dproba(indwrong(j)),'%1.3f'),' p-value = ', num2str(proba(indwrong(j)),'%1.3f')])
end
skipline()
+ indcorr1 = stab_map_2(lpmat(iwrong,:),alpha2, pvalue_corr, awrongname);
+ indcorr2 = stab_map_2(lpmat(itmp,:),alpha2, pvalue_corr, ['NOT_',awrongname]);
+ indcorr = union(indcorr1(:), indcorr2(:));
+ indcorr = indcorr(~ismember(indcorr(:),indwrong));
+ indwrong = [indwrong(:); indcorr(:)];
+ skipline()
if ~isempty(indwrong)
- stab_map_1(lpmat, [1:Nsam], iwrong, awronguniname, 1, indwrong, OutputDirectoryName,[],awrongunititle);
+% stab_map_1(lpmat, itmp, iwrong, awronguniname, 1, indwrong, OutputDirectoryName,[],awrongunititle);
+ scatter_mcf(lpmat(itmp,indunst),lpmat(ixun,indunst),cellstr(bayestopt_.name(indunst+nshock)), ...
+ '.', [fname_,'_',awrongname], OutputDirectoryName,awrongtitle,[], options_, ...
+ 'NO inability to fund a solution', 'inability to find a solution')
end
end
- skipline()
- disp('Starting bivariate analysis:')
-
- c0=corrcoef(lpmat(istable,:));
- c00=tril(c0,-1);
-
- 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
- if length(iindeterm)>10,
- stab_map_2(lpmat(iindeterm,:),alpha2, pvalue_corr, aindname, OutputDirectoryName,xparam1,aindtitle);
- end
- if length(ixun)>10,
- stab_map_2(lpmat(ixun,:),alpha2, pvalue_corr, aunstname, OutputDirectoryName,xparam1,aunsttitle);
- end
- if length(iwrong)>10,
- stab_map_2(lpmat(iwrong,:),alpha2, pvalue_corr, awrongname, OutputDirectoryName,xparam1,awrongtitle);
+ if ~isempty(irestriction),
+ [proba, dproba] = stab_map_1(lpmat, irestriction, inorestriction, awronguniname,0);
+ % indwrong=find(dproba>ksstat);
+ indcalib=find(proba<pvalue_ks);
+ [tmp,jtmp] = sort(proba(indcalib),2,'ascend');
+ indcalib = indcalib(jtmp);
+ disp('Smirnov statistics in driving prior restrictions')
+ for j=1:length(indcalib),
+ disp([M_.param_names(estim_params_.param_vals(indcalib(j),1),:),' d-stat = ', num2str(dproba(indcalib(j)),'%1.3f'),' p-value = ', num2str(proba(indcalib(j)),'%1.3f')])
+ end
+ skipline()
+ indcorr1 = stab_map_2(lpmat(irestriction,:),alpha2, pvalue_corr, acalibname);
+ indcorr2 = stab_map_2(lpmat(inorestriction,:),alpha2, pvalue_corr, ['NOT_',acalibname]);
+ indcorr = union(indcorr1(:), indcorr2(:));
+ indcorr = indcorr(~ismember(indcorr(:),indcalib));
+ indcalib = [indcalib(:); indcorr(:)];
+ skipline()
+ if ~isempty(indcalib)
+% stab_map_1(lpmat, irestriction, inorestriction, acalibuniname, 1, indcalib, OutputDirectoryName,[],acalibunititle);
+ scatter_mcf(lpmat(irestriction,indcalib),lpmat(inorestriction,indcalib),cellstr(bayestopt_.name(indcalib+nshock)), ...
+ '.', [fname_,'_',acalibname], OutputDirectoryName,acalibtitle,[], options_, ...
+ 'prior IRF/moment calibration', 'NO prior IRF/moment calibration')
+ end
+ iok = irestriction(1);
+ else
+ iok = istable(1);
end
+% skipline()
+% disp('Starting bivariate analysis:')
+%
+% c0=corrcoef(lpmat(istable,:));
+% c00=tril(c0,-1);
+%
+% 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
+% if length(iindeterm)>10,
+% stab_map_2(lpmat(iindeterm,:),alpha2, pvalue_corr, aindname, OutputDirectoryName,xparam1,aindtitle);
+% end
+% if length(ixun)>10,
+% stab_map_2(lpmat(ixun,:),alpha2, pvalue_corr, aunstname, OutputDirectoryName,xparam1,aunsttitle);
+% end
+% if length(iwrong)>10,
+% stab_map_2(lpmat(iwrong,:),alpha2, pvalue_corr, awrongname, OutputDirectoryName,xparam1,awrongtitle);
+% end
+% if length(irestriction)>10,
+% stab_map_2(lpmat(irestriction,:),alpha2, pvalue_corr, acalibname, OutputDirectoryName,xparam1,acalibtitle);
+% iok = irestriction(1);
+% else
+% iok = istable(1);
+% end
x0=0.5.*(bounds.ub(1:nshock)-bounds.lb(1:nshock))+bounds.lb(1:nshock);
- x0 = [x0; lpmat(istable(1),:)'];
- if istable(end)~=Nsam
- M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(1),:)';
- [oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
- % stoch_simul([]);
- end
+ x0 = [x0; lpmat(iok,:)'];
+ M_.params(estim_params_.param_vals(:,1)) = lpmat(iok,:)';
+ [oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
+ % stoch_simul([]);
else
disp('All parameter values in the specified ranges are not acceptable!')
x0=[];
end
else
- disp('All parameter values in the specified ranges give unique saddle-path solution!')
- x0=0.5.*(bounds.ub(1:nshock)-bounds.lb(1:nshock))+bounds.lb(1:nshock);
- x0 = [x0; lpmat(istable(1),:)'];
+ disp('All parameter values in the specified ranges give unique saddle-path solution,')
+ disp('and match prior IRF/moment restriction(s) if any!')
+ x0=0.5.*(bounds.ub(1:nshock)-bounds.lb(1:nshock))+bounds.lb(1:nshock);
+ x0 = [x0; lpmat(istable(1),:)'];
end
diff --git a/matlab/gsa/stab_map_2.m b/matlab/gsa/stab_map_2.m
index d4e94f02e82221e4a5d4338f1dfb213e0af31a31..7848de1c6e235bb3559aa82342693ba1a10aa174 100644
--- a/matlab/gsa/stab_map_2.m
+++ b/matlab/gsa/stab_map_2.m
@@ -1,4 +1,4 @@
-function stab_map_2(x,alpha2, pvalue, fnam, dirname,xparam1,figtitle)
+function indcorr = stab_map_2(x,alpha2, pvalue_crit, fnam, dirname,xparam1,figtitle)
% function stab_map_2(x, alpha2, pvalue, fnam, dirname,xparam1)
%
% Written by Marco Ratto
@@ -32,11 +32,13 @@ global bayestopt_ estim_params_ options_ oo_ M_
npar=size(x,2);
nsam=size(x,1);
ishock= npar>estim_params_.np;
+nograph = options_.nograph;
if nargin<4,
fnam='';
end
if nargin<5,
dirname='';
+ nograph=1;
end
if nargin<6,
xparam1=[];
@@ -53,7 +55,7 @@ nshock = nshock + estim_params_.nvn;
nshock = nshock + estim_params_.ncx;
nshock = nshock + estim_params_.ncn;
-c0=corrcoef(x);
+[c0, pvalue] = corrcoef(x);
c00=tril(c0,-1);
fig_nam_=[fname_,'_',fnam,'_corr_'];
@@ -70,13 +72,13 @@ end
disp([' '])
disp(['Correlation analysis for ',fnam])
+indcorr = [];
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,
+ if pvalue(j,i2(jx))<pvalue_crit,
+ indcorr = [indcorr; [j i2(jx)]];
j2=j2+1;
if ishock,
tmp_name = (['[',bayestopt_.name{j},',',bayestopt_.name{i2(jx)},']']);
@@ -85,7 +87,7 @@ for j=1:npar,
end
fprintf(1,'%20s: corrcoef = %7.3f\n',tmp_name,c0(i2(jx),j));
- if ~options_.nograph,
+ if ~nograph,
if mod(j2,12)==1,
ifig=ifig+1;
hh=dyn_figure(options_,'name',[figtitle,' sample bivariate projection ', num2str(ifig)]);
@@ -118,12 +120,12 @@ for j=1:npar,
end
end
- if ~options_.nograph && (j==(npar)) && j2>0 && (mod(j2,12)~=0),
+ if ~nograph && (j==(npar)) && j2>0 && (mod(j2,12)~=0),
dyn_saveas(hh,[dirname,filesep,fig_nam_,int2str(ifig)],options_);
end
end
if j2==0,
- disp(['No correlation term with pvalue <', num2str(pvalue),' found for ',fnam])
+ disp(['No correlation term with pvalue <', num2str(pvalue_crit),' and |corr. coef.| >',num2str(alpha2),' found for ',fnam])
end
%close all