stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Identification: various cosmetic fixes

new-samplers
Johannes Pfeifer 2023-11-30 12:26:58 +01:00 committed by Sébastien Villemot
parent e68793030c
commit 392721097c
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
3 changed files with 58 additions and 51 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
matlab/dynare_identification.m
View File

@ -1,5 +1,5 @@
function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = dynare_identification(M_,oo_,options_,bayestopt_,estim_params_,options_ident, pdraws0) function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = dynare_identification(M_,oo_,options_,bayestopt_,estim_params_,options_ident, pdraws0)
%function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = dynare_identification(options_ident, pdraws0) % [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = dynare_identification(options_ident, pdraws0)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called, when the user specifies identification(...); in the mod file. It prepares all identification analysis: % This function is called, when the user specifies identification(...); in the mod file. It prepares all identification analysis:
% (1) set options, local and persistent variables for a new identification % (1) set options, local and persistent variables for a new identification
@ -11,6 +11,11 @@ function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, ST
% to put identification in your mod file, otherwise the preprocessor won't provide all necessary objects % to put identification in your mod file, otherwise the preprocessor won't provide all necessary objects
% ========================================================================= % =========================================================================
% INPUTS % INPUTS
% * M_ [structure] Matlab's structure describing the model
% * oo_ [structure] Matlab's structure describing the results
% * options_ [structure] Matlab's structure describing the current options
% * bayestopt_ [structure] describing the priors
% * estim_params_ [structure] characterizing parameters to be estimated
% * options_ident [structure] identification options % * options_ident [structure] identification options
% * pdraws0 [SampleSize by totparam_nbr] optional: matrix of MC sample of model parameters % * pdraws0 [SampleSize by totparam_nbr] optional: matrix of MC sample of model parameters
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
@ -781,6 +786,19 @@ if iload <=0
else else
maxrun_dMINIMAL = 0; maxrun_dMINIMAL = 0;
end end
si_dDYNAMICnorm=NaN(max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL]),size(STO_si_dDYNAMIC,2));
if ~options_MC.no_identification_reducedform
si_dREDUCEDFORMnorm=NaN(max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL]),size(STO_si_dREDUCEDFORM,2));
end
if ~options_MC.no_identification_moments
si_dMOMENTSnorm=NaN(max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL]),size(STO_si_dMOMENTS,2));
end
if ~options_MC.no_identification_spectrum
dSPECTRUMnorm=NaN(max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL]),size(STO_dSPECTRUM,2));
end
if ~options_MC.no_identification_minimal
dMINIMALnorm=NaN(max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL]),size(STO_dMINIMAL,2));
end
for irun=1:max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL]) for irun=1:max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL])
iter=iter+1; iter=iter+1;
% note that this is not the same si_dDYNAMICnorm as computed in identification_analysis % note that this is not the same si_dDYNAMICnorm as computed in identification_analysis
@ -912,7 +930,7 @@ if SampleSize > 1
fprintf('Testing %s.\n',tittxt); fprintf('Testing %s.\n',tittxt);
if ~iload if ~iload
options_ident.tittxt = tittxt; %title text for graphs and figures options_ident.tittxt = tittxt; %title text for graphs and figures
[ide_moments_min, ide_spectrum_min, ide_minimal_min, ide_hess_min, ide_reducedform_min, ide_dynamic_min, derivatives_info_min, info_min, error_indicator_min] = ... [ide_moments_min, ide_spectrum_min, ide_minimal_min, ide_hess_min, ide_reducedform_min, ide_dynamic_min, ~, ~, error_indicator_min] = ...
identification_analysis(M_,options_,oo_,bayestopt_,estim_params_,pdraws(jmin,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes persistent variables identification_analysis(M_,options_,oo_,bayestopt_,estim_params_,pdraws(jmin,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes persistent variables
save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_min', 'ide_moments_min','ide_spectrum_min','ide_minimal_min','ide_reducedform_min', 'ide_dynamic_min', 'jmin', '-append'); save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_min', 'ide_moments_min','ide_spectrum_min','ide_minimal_min','ide_reducedform_min', 'ide_dynamic_min', 'jmin', '-append');
end end

4
matlab/identification_analysis.m
View File

@ -142,7 +142,7 @@ error_indicator.identification_spectrum=0;
if info(1) == 0 %no errors in solution if info(1) == 0 %no errors in solution
% Compute parameter Jacobians for identification analysis % Compute parameter Jacobians for identification analysis
[MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params_, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state); [~, ~, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params_, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if isempty(dMINIMAL) if isempty(dMINIMAL)
% Komunjer and Ng is not computed if (1) minimality conditions are not fullfilled or (2) there are more shocks and measurement errors than observables, so we need to reset options % Komunjer and Ng is not computed if (1) minimality conditions are not fullfilled or (2) there are more shocks and measurement errors than observables, so we need to reset options
error_indicator.identification_minimal = 1; error_indicator.identification_minimal = 1;
@ -300,7 +300,7 @@ if info(1) == 0 %no errors in solution
derivatives_info.no_DLIK = 1; derivatives_info.no_DLIK = 1;
bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding
%note that for order>1 we do not provide any information on DT,DYss,DOM in derivatives_info, such that dsge_likelihood creates an error. Therefore the computation will be based on simulated_moment_uncertainty for order>1. %note that for order>1 we do not provide any information on DT,DYss,DOM in derivatives_info, such that dsge_likelihood creates an error. Therefore the computation will be based on simulated_moment_uncertainty for order>1.
[fval, info, cost_flag, DLIK, AHess, ys, trend_coeff, M_, options_, bayestopt_, oo_.dr] = dsge_likelihood(params', dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr, oo_.steady_state,oo_.exo_steady_state, oo_.exo_det_steady_state. derivatives_info); %non-used output variables need to be set for octave for some reason [~, info, ~, ~, AHess, ~, ~, M_, options_, ~, oo_.dr] = dsge_likelihood(params', dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr, oo_.steady_state,oo_.exo_steady_state, oo_.exo_det_steady_state. derivatives_info); %non-used output variables need to be set for octave for some reason
%note that for the order of parameters in AHess we have: stderr parameters come first, corr parameters second, model parameters third. the order within these blocks corresponds to the order specified in the estimated_params block %note that for the order of parameters in AHess we have: stderr parameters come first, corr parameters second, model parameters third. the order within these blocks corresponds to the order specified in the estimated_params block
options_.analytic_derivation = analytic_derivation; %reset option options_.analytic_derivation = analytic_derivation; %reset option
AHess = -AHess; %take negative of hessian AHess = -AHess; %take negative of hessian

83
matlab/plot_identification.m
View File

@ -1,5 +1,5 @@
function plot_identification(M_, params, idemoments, idehess, idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, fname, options_, estim_params_, bayestopt_, tit_TeX, name_tex) function plot_identification(M_, params, idemoments, idehess, idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, fname, options_, estim_params_, bayestopt_, tit_TeX, name_tex)
% function plot_identification(params,idemoments,idehess,idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, fname, options_, estim_params_, bayestopt_, tit_TeX, name_tex) % plot_identification(M_, params,idemoments,idehess,idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, fname, options_, estim_params_, bayestopt_, tit_TeX, name_tex)
% %
% INPUTS % INPUTS
% o M_ [structure] model % o M_ [structure] model
@ -10,10 +10,15 @@ function plot_identification(M_, params, idemoments, idehess, idemodel, idelre,
% o idelre [structure] identification results for the LRE model % o idelre [structure] identification results for the LRE model
% o advanced [integer] flag for advanced identification checks % o advanced [integer] flag for advanced identification checks
% o tittxt [char] name of the results to plot % o tittxt [char] name of the results to plot
% o name [char] list of names % o name [char] list of parameter names
% o IdentifDirectoryName [char] directory name % o IdentifDirectoryName [char] directory name
% o fname [char] file name
% o options_ [structure] structure describing the current options
% o estim_params_ [structure] characterizing parameters to be estimated
% o bayestopt_ [structure] describing the priors
% o tittxt [char] TeX-name of the results to plot % o tittxt [char] TeX-name of the results to plot
% o name_tex [char] TeX-names of the parameters % o name_tex [char] TeX-names of the parameters
%
% OUTPUTS % OUTPUTS
% None % None
% %
@ -53,20 +58,19 @@ si_dLREnorm = idelre.si_dDYNAMICnorm;
tittxt1=regexprep(tittxt, ' ', '_'); tittxt1=regexprep(tittxt, ' ', '_');
tittxt1=strrep(tittxt1, '.', ''); tittxt1=strrep(tittxt1, '.', '');
if SampleSize == 1 if SampleSize == 1
si_dMOMENTS = idemoments.si_dMOMENTS;
hh_fig = dyn_figure(options_.nodisplay,'Name',[tittxt, ' - Identification using info from observables']); hh_fig = dyn_figure(options_.nodisplay,'Name',[tittxt, ' - Identification using info from observables']);
subplot(211) subplot(211)
mmm = (idehess.ide_strength_dMOMENTS); mmm = (idehess.ide_strength_dMOMENTS);
[ss, is] = sort(mmm); [~, is] = sort(mmm);
if ~all(isnan(idehess.ide_strength_dMOMENTS_prior)) ... if ~all(isnan(idehess.ide_strength_dMOMENTS_prior)) ...
&& ~(nparam == 1 && ~isoctave && matlab_ver_less_than('9.7')) % MATLAB < R2019b does not accept bar(1, [2 3]) && ~(nparam == 1 && ~isoctave && matlab_ver_less_than('9.7')) % MATLAB < R2019b does not accept bar(1, [2 3])
bar(1:nparam,log([idehess.ide_strength_dMOMENTS(:,is)' idehess.ide_strength_dMOMENTS_prior(:,is)'])) bar(1:nparam,log([idehess.ide_strength_dMOMENTS(:,is)' idehess.ide_strength_dMOMENTS_prior(:,is)']))
else else
bar(1:nparam,log([idehess.ide_strength_dMOMENTS(:,is)' ])) bar(1:nparam,log(idehess.ide_strength_dMOMENTS(:,is)'))
end end
hold on hold on
plot((1:length(idehess.ide_strength_dMOMENTS(:,is)))-0.15,log([idehess.ide_strength_dMOMENTS(:,is)']),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none') plot((1:length(idehess.ide_strength_dMOMENTS(:,is)))-0.15,log(idehess.ide_strength_dMOMENTS(:,is)'),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none')
plot((1:length(idehess.ide_strength_dMOMENTS_prior(:,is)))+0.15,log([idehess.ide_strength_dMOMENTS_prior(:,is)']),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none') plot((1:length(idehess.ide_strength_dMOMENTS_prior(:,is)))+0.15,log(idehess.ide_strength_dMOMENTS_prior(:,is)'),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none')
if any(isinf(log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices)))) if any(isinf(log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices))))
%-Inf, i.e. 0 strength %-Inf, i.e. 0 strength
inf_indices=find(isinf(log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices))) & log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices))<0); inf_indices=find(isinf(log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices))) & log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices))<0);
@ -94,8 +98,8 @@ if SampleSize == 1
if options_.TeX if options_.TeX
text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex') text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else else
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end end
end end
if ~all(isnan(idehess.ide_strength_dMOMENTS_prior)) if ~all(isnan(idehess.ide_strength_dMOMENTS_prior))
legend('relative to param value','relative to prior std','Location','Best') legend('relative to param value','relative to prior std','Location','Best')
@ -116,8 +120,8 @@ if SampleSize == 1
bar(1:nparam, log([idehess.deltaM(is)])) bar(1:nparam, log([idehess.deltaM(is)]))
end end
hold on hold on
plot((1:length(idehess.deltaM(is)))-0.15,log([idehess.deltaM(is)']),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none') plot((1:length(idehess.deltaM(is)))-0.15,log(idehess.deltaM(is)'),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none')
plot((1:length(idehess.deltaM_prior(is)))+0.15,log([idehess.deltaM_prior(is)']),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none') plot((1:length(idehess.deltaM_prior(is)))+0.15,log(idehess.deltaM_prior(is)'),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none')
inf_pos=find(isinf(log(idehess.deltaM))); inf_pos=find(isinf(log(idehess.deltaM)));
if ~isempty(inf_pos) if ~isempty(inf_pos)
inf_indices=~ismember(inf_pos,idehess.sensitivity_zero_pos); inf_indices=~ismember(inf_pos,idehess.sensitivity_zero_pos);
@ -137,8 +141,8 @@ if SampleSize == 1
if options_.TeX if options_.TeX
text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex') text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else else
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end end
end end
if ~all(isnan(idehess.deltaM_prior)) if ~all(isnan(idehess.deltaM_prior))
legend('relative to param value','relative to prior std','Location','Best') legend('relative to param value','relative to prior std','Location','Best')
@ -191,8 +195,8 @@ if SampleSize == 1
if options_.TeX if options_.TeX
text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex') text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else else
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end end
end end
legend('Moments','Model','LRE model','Location','Best') legend('Moments','Model','LRE model','Location','Best')
title('Sensitivity bars using derivatives (log-scale)') title('Sensitivity bars using derivatives (log-scale)')
@ -214,9 +218,6 @@ if SampleSize == 1
% identificaton patterns % identificaton patterns
for j=1:size(idemoments.cosndMOMENTS,2) for j=1:size(idemoments.cosndMOMENTS,2)
pax=NaN(nparam,nparam); pax=NaN(nparam,nparam);
% fprintf('\n')
% disp(['Collinearity patterns with ', int2str(j) ,' parameter(s)'])
% fprintf('%-15s [%-*s] %10s\n','Parameter',(15+1)*j,' Expl. params ','cosn')
for i=1:nparam for i=1:nparam
namx=''; namx='';
for in=1:j for in=1:j
@ -227,12 +228,11 @@ if SampleSize == 1
if options_.TeX if options_.TeX
namx=[namx ' ' sprintf('%-15s',name_tex{dumpindx})]; namx=[namx ' ' sprintf('%-15s',name_tex{dumpindx})];
else else
namx=[namx ' ' sprintf('%-15s',name{dumpindx})]; namx=[namx ' ' sprintf('%-15s',name{dumpindx})];
end end
pax(i,dumpindx)=idemoments.cosndMOMENTS(i,j); pax(i,dumpindx)=idemoments.cosndMOMENTS(i,j);
end end
end end
% fprintf('%-15s [%s] %10.3f\n',name{i},namx,idemoments.cosndMOMENTS(i,j))
end end
hh_fig = dyn_figure(options_.nodisplay,'Name',[tittxt,' - Collinearity patterns with ', int2str(j) ,' parameter(s)']); hh_fig = dyn_figure(options_.nodisplay,'Name',[tittxt,' - Collinearity patterns with ', int2str(j) ,' parameter(s)']);
imagesc(pax,[0 1]); imagesc(pax,[0 1]);
@ -243,9 +243,9 @@ if SampleSize == 1
text(ip,(0.5),name_tex{ip},'rotation',90,'HorizontalAlignment','left','interpreter','latex') text(ip,(0.5),name_tex{ip},'rotation',90,'HorizontalAlignment','left','interpreter','latex')
text(0.5,ip,name_tex{ip},'rotation',0,'HorizontalAlignment','right','interpreter','latex') text(0.5,ip,name_tex{ip},'rotation',0,'HorizontalAlignment','right','interpreter','latex')
else else
text(ip,(0.5),name{ip},'rotation',90,'HorizontalAlignment','left','interpreter','none') text(ip,(0.5),name{ip},'rotation',90,'HorizontalAlignment','left','interpreter','none')
text(0.5,ip,name{ip},'rotation',0,'HorizontalAlignment','right','interpreter','none') text(0.5,ip,name{ip},'rotation',0,'HorizontalAlignment','right','interpreter','none')
end end
end end
colorbar; colorbar;
colormap('jet'); colormap('jet');
@ -275,7 +275,7 @@ if SampleSize == 1
end end
end end
skipline() skipline()
[U,S,V]=svd(idehess.AHess,0); [~,S,V]=svd(idehess.AHess,0);
S=diag(S); S=diag(S);
if idehess.flag_score if idehess.flag_score
if nparam<5 if nparam<5
@ -288,8 +288,6 @@ if SampleSize == 1
tex_tit_2=[tittxt,' - Identification patterns (Information matrix): HIGHEST SV']; tex_tit_2=[tittxt,' - Identification patterns (Information matrix): HIGHEST SV'];
end end
else else
% S = idemoments.S;
% V = idemoments.V;
if nparam<5 if nparam<5
f1 = dyn_figure(options_.nodisplay,'Name',[tittxt,' - Identification patterns (moments Information matrix)']); f1 = dyn_figure(options_.nodisplay,'Name',[tittxt,' - Identification patterns (moments Information matrix)']);
tex_tit_1=[tittxt,' - Identification patterns (moments Information matrix)']; tex_tit_1=[tittxt,' - Identification patterns (moments Information matrix)'];
@ -322,10 +320,10 @@ if SampleSize == 1
if options_.TeX if options_.TeX
text(ip,-0.02,name_tex{ip},'rotation',90,'HorizontalAlignment','right','interpreter','latex') text(ip,-0.02,name_tex{ip},'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else else
text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none') text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end end
end end
end
title(['Singular value ',num2str(Stit)]) title(['Singular value ',num2str(Stit)])
end end
dyn_saveas(f1,[ IdentifDirectoryName '/' fname '_ident_pattern_' tittxt1 '_1' ],options_.nodisplay,options_.graph_format); dyn_saveas(f1,[ IdentifDirectoryName '/' fname '_ident_pattern_' tittxt1 '_1' ],options_.nodisplay,options_.graph_format);
@ -361,10 +359,10 @@ if SampleSize == 1
end end
else else
hh_fig = dyn_figure(options_.nodisplay,'Name',['MC sensitivities']); hh_fig = dyn_figure(options_.nodisplay,'Name','MC sensitivities');
subplot(211) subplot(211)
mmm = (idehess.ide_strength_dMOMENTS); mmm = (idehess.ide_strength_dMOMENTS);
[ss, is] = sort(mmm); [~, is] = sort(mmm);
mmm = mean(si_dMOMENTSnorm)'; mmm = mean(si_dMOMENTSnorm)';
mmm = mmm./max(mmm); mmm = mmm./max(mmm);
if advanced if advanced
@ -384,8 +382,8 @@ else
if options_.TeX if options_.TeX
text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex') text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else else
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end end
end end
if advanced if advanced
legend('Moments','Model','LRE model','Location','Best') legend('Moments','Model','LRE model','Location','Best')
@ -407,11 +405,10 @@ else
end end
if advanced if advanced
if ~options_.nodisplay, if ~options_.nodisplay
skipline() skipline()
disp('Displaying advanced diagnostics') disp('Displaying advanced diagnostics')
end end
% options_.nograph=1;
hh_fig = dyn_figure(options_.nodisplay,'Name','MC Condition Number'); hh_fig = dyn_figure(options_.nodisplay,'Name','MC Condition Number');
subplot(221) subplot(221)
hist(log10(idemodel.cond)) hist(log10(idemodel.cond))
@ -448,16 +445,6 @@ else
options_mcf.title = 'MC Highest Condition Number Model Moments'; options_mcf.title = 'MC Highest Condition Number Model Moments';
[~,is]=sort(idemoments.cond); [~,is]=sort(idemoments.cond);
mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_); mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
% [proba, dproba] = stab_map_1(idemoments.Mco', is(1:ncut), is(ncut+1:end), 'HighestCondNumberMoments_vs_Mco', 1, [], IdentifDirectoryName);
% for j=1:nparam,
% % ibeh=find(idemoments.Mco(j,:)<0.9);
% % inonbeh=find(idemoments.Mco(j,:)>=0.9);
% % if ~isempty(ibeh) && ~isempty(inonbeh)
% % [proba, dproba] = stab_map_1(params, ibeh, inonbeh, ['HighestMultiCollinearity_',name{j}], 1, [], IdentifDirectoryName);
% % end
% [~,is]=sort(idemoments.Mco(:,j));
% [proba, dproba] = stab_map_1(params, is(1:ncut), is(ncut+1:end), ['MC_HighestMultiCollinearity_',name{j}], 1, [], IdentifDirectoryName, 0.15);
% end
if nparam<5 if nparam<5
f1 = dyn_figure(options_.nodisplay,'Name',[tittxt,' - MC Identification patterns (moments): HIGHEST SV']); f1 = dyn_figure(options_.nodisplay,'Name',[tittxt,' - MC Identification patterns (moments): HIGHEST SV']);
@ -487,8 +474,10 @@ else
SSS = idemoments.S(:,jj); SSS = idemoments.S(:,jj);
end end
subplot(nsubplo,1,jj) subplot(nsubplo,1,jj)
post_median=NaN(1,nparam);
hpd_interval=NaN(nparam,2);
for i=1:nparam for i=1:nparam
[post_mean, post_median(:,i), post_var, hpd_interval(i,:), post_deciles] = posterior_moments(VVV(:,i),0,0.9); [~, post_median(:,i), ~, hpd_interval(i,:)] = posterior_moments(VVV(:,i),0,0.9);
end end
bar(post_median) bar(post_median)
hold on, plot(hpd_interval,'--*r'), hold on, plot(hpd_interval,'--*r'),
@ -500,10 +489,10 @@ else
if options_.TeX if options_.TeX
text(ip,-0.02,name_tex{ip},'rotation',90,'HorizontalAlignment','right','interpreter','latex') text(ip,-0.02,name_tex{ip},'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else else
text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none') text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end end
end end
end
title(['MEAN Singular value ',num2str(Stit)]) title(['MEAN Singular value ',num2str(Stit)])
end end
dyn_saveas(f1,[IdentifDirectoryName '/' fname '_MC_ident_pattern_1' ],options_.nodisplay,options_.graph_format); dyn_saveas(f1,[IdentifDirectoryName '/' fname '_MC_ident_pattern_1' ],options_.nodisplay,options_.graph_format);