diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m index b733bb74d..2af64d33f 100644 --- a/matlab/dynare_identification.m +++ b/matlab/dynare_identification.m @@ -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(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: % (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 % ========================================================================= % 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 % * pdraws0 [SampleSize by totparam_nbr] optional: matrix of MC sample of model parameters % ------------------------------------------------------------------------- @@ -781,6 +786,19 @@ if iload <=0 else maxrun_dMINIMAL = 0; 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]) iter=iter+1; % 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); if ~iload 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 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 diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m index 42ddaffe3..09c9cefe3 100644 --- a/matlab/identification_analysis.m +++ b/matlab/identification_analysis.m @@ -142,7 +142,7 @@ error_indicator.identification_spectrum=0; if info(1) == 0 %no errors in solution % 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) % 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; @@ -300,7 +300,7 @@ if info(1) == 0 %no errors in solution 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 %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 options_.analytic_derivation = analytic_derivation; %reset option AHess = -AHess; %take negative of hessian diff --git a/matlab/plot_identification.m b/matlab/plot_identification.m index fcfbc6074..19b7836ae 100644 --- a/matlab/plot_identification.m +++ b/matlab/plot_identification.m @@ -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(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 % 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 advanced [integer] flag for advanced identification checks % 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 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 name_tex [char] TeX-names of the parameters +% % OUTPUTS % None % @@ -53,20 +58,19 @@ si_dLREnorm = idelre.si_dDYNAMICnorm; tittxt1=regexprep(tittxt, ' ', '_'); tittxt1=strrep(tittxt1, '.', ''); if SampleSize == 1 - si_dMOMENTS = idemoments.si_dMOMENTS; hh_fig = dyn_figure(options_.nodisplay,'Name',[tittxt, ' - Identification using info from observables']); subplot(211) mmm = (idehess.ide_strength_dMOMENTS); - [ss, is] = sort(mmm); + [~, is] = sort(mmm); 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]) bar(1:nparam,log([idehess.ide_strength_dMOMENTS(:,is)' idehess.ide_strength_dMOMENTS_prior(:,is)'])) else - bar(1:nparam,log([idehess.ide_strength_dMOMENTS(:,is)' ])) + bar(1:nparam,log(idehess.ide_strength_dMOMENTS(:,is)')) end 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_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(:,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') if any(isinf(log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices)))) %-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); @@ -94,8 +98,8 @@ if SampleSize == 1 if options_.TeX text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex') else - text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') - end + text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') + end end if ~all(isnan(idehess.ide_strength_dMOMENTS_prior)) 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)])) end 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_prior(is)))+0.15,log([idehess.deltaM_prior(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') inf_pos=find(isinf(log(idehess.deltaM))); if ~isempty(inf_pos) inf_indices=~ismember(inf_pos,idehess.sensitivity_zero_pos); @@ -137,8 +141,8 @@ if SampleSize == 1 if options_.TeX text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex') else - text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') - end + text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') + end end if ~all(isnan(idehess.deltaM_prior)) legend('relative to param value','relative to prior std','Location','Best') @@ -191,8 +195,8 @@ if SampleSize == 1 if options_.TeX text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex') else - text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') - end + text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') + end end legend('Moments','Model','LRE model','Location','Best') title('Sensitivity bars using derivatives (log-scale)') @@ -214,9 +218,6 @@ if SampleSize == 1 % identificaton patterns for j=1:size(idemoments.cosndMOMENTS,2) 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 namx=''; for in=1:j @@ -227,12 +228,11 @@ if SampleSize == 1 if options_.TeX namx=[namx ' ' sprintf('%-15s',name_tex{dumpindx})]; else - namx=[namx ' ' sprintf('%-15s',name{dumpindx})]; + namx=[namx ' ' sprintf('%-15s',name{dumpindx})]; end pax(i,dumpindx)=idemoments.cosndMOMENTS(i,j); end end - % fprintf('%-15s [%s] %10.3f\n',name{i},namx,idemoments.cosndMOMENTS(i,j)) end hh_fig = dyn_figure(options_.nodisplay,'Name',[tittxt,' - Collinearity patterns with ', int2str(j) ,' parameter(s)']); 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(0.5,ip,name_tex{ip},'rotation',0,'HorizontalAlignment','right','interpreter','latex') else - text(ip,(0.5),name{ip},'rotation',90,'HorizontalAlignment','left','interpreter','none') - text(0.5,ip,name{ip},'rotation',0,'HorizontalAlignment','right','interpreter','none') - end + text(ip,(0.5),name{ip},'rotation',90,'HorizontalAlignment','left','interpreter','none') + text(0.5,ip,name{ip},'rotation',0,'HorizontalAlignment','right','interpreter','none') + end end colorbar; colormap('jet'); @@ -275,7 +275,7 @@ if SampleSize == 1 end end skipline() - [U,S,V]=svd(idehess.AHess,0); + [~,S,V]=svd(idehess.AHess,0); S=diag(S); if idehess.flag_score if nparam<5 @@ -288,8 +288,6 @@ if SampleSize == 1 tex_tit_2=[tittxt,' - Identification patterns (Information matrix): HIGHEST SV']; end else - % S = idemoments.S; - % V = idemoments.V; if nparam<5 f1 = dyn_figure(options_.nodisplay,'Name',[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 text(ip,-0.02,name_tex{ip},'rotation',90,'HorizontalAlignment','right','interpreter','latex') 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 title(['Singular value ',num2str(Stit)]) end dyn_saveas(f1,[ IdentifDirectoryName '/' fname '_ident_pattern_' tittxt1 '_1' ],options_.nodisplay,options_.graph_format); @@ -361,10 +359,10 @@ if SampleSize == 1 end else - hh_fig = dyn_figure(options_.nodisplay,'Name',['MC sensitivities']); + hh_fig = dyn_figure(options_.nodisplay,'Name','MC sensitivities'); subplot(211) mmm = (idehess.ide_strength_dMOMENTS); - [ss, is] = sort(mmm); + [~, is] = sort(mmm); mmm = mean(si_dMOMENTSnorm)'; mmm = mmm./max(mmm); if advanced @@ -384,8 +382,8 @@ else if options_.TeX text(ip,dy(1),name_tex{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','latex') else - text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') - end + text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') + end end if advanced legend('Moments','Model','LRE model','Location','Best') @@ -407,11 +405,10 @@ else end if advanced - if ~options_.nodisplay, + if ~options_.nodisplay skipline() disp('Displaying advanced diagnostics') end - % options_.nograph=1; hh_fig = dyn_figure(options_.nodisplay,'Name','MC Condition Number'); subplot(221) hist(log10(idemodel.cond)) @@ -448,16 +445,6 @@ else options_mcf.title = 'MC Highest Condition Number Model Moments'; [~,is]=sort(idemoments.cond); 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 f1 = dyn_figure(options_.nodisplay,'Name',[tittxt,' - MC Identification patterns (moments): HIGHEST SV']); @@ -487,8 +474,10 @@ else SSS = idemoments.S(:,jj); end subplot(nsubplo,1,jj) + post_median=NaN(1,nparam); + hpd_interval=NaN(nparam,2); 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 bar(post_median) hold on, plot(hpd_interval,'--*r'), @@ -500,10 +489,10 @@ else if options_.TeX text(ip,-0.02,name_tex{ip},'rotation',90,'HorizontalAlignment','right','interpreter','latex') 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 title(['MEAN Singular value ',num2str(Stit)]) end dyn_saveas(f1,[IdentifDirectoryName '/' fname '_MC_ident_pattern_1' ],options_.nodisplay,options_.graph_format);