diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m index 9ec45a30e..fb0ed930e 100644 --- a/matlab/dynare_identification.m +++ b/matlab/dynare_identification.m @@ -78,7 +78,7 @@ options_.Schur_vec_tol = 1.e-8; options_ = set_default_option(options_,'datafile',[]); options_.mode_compute = 0; options_.plot_priors = 0; -[data,rawdata,xparam1,data_info]=dynare_estimation_init([],fname_,1); +[data,rawdata,xparam1,data_info]=dynare_estimation_init(M_.endo_names,fname_,1); if isempty(data_info), data_info.gend = periods; data_info.data = []; @@ -138,6 +138,10 @@ else name_tex = [cellstr(M_.exo_names_tex); cellstr(M_.param_names_tex)]; end +options_ident = set_default_option(options_ident,'max_ord_bruteforce',min([2,nparam-1])); +max_ord_bruteforce = options_ident.max_ord_bruteforce; + + MaxNumberOfBytes=options_.MaxNumberOfBytes; disp(' ') @@ -324,12 +328,14 @@ if iload <=0, identification_checks(H(indH,:)./normH(:,ones(nparam,1)),JJ(indJJ,:)./normJ(:,ones(nparam,1)), gp(indLRE,:)./normLRE(:,ones(size(gp,2),1))); % identification_checks(H(indH,:),JJ(indJJ,:), gp(indLRE,:), bayestopt_); indok = find(max(idemoments.indno{iteration},[],1)==0); - if iteration ==1 && ~isempty(indok), + if iteration ==1, ide_strength_J=NaN(1,nparam); ide_strength_J_prior=NaN(1,nparam); - if advanced, - [pars, cosnJ] = ident_bruteforce(JJ(indJJ,:)./normJ(:,ones(nparam,1)),2,1,name_tex); - end + end + if iteration ==1 && advanced, + [pars, cosnJ] = ident_bruteforce(JJ(indJJ,:)./normJ(:,ones(nparam,1)),max_ord_bruteforce,options_.TeX,name_tex); + end + if iteration ==1 && ~isempty(indok), normaliz = abs(params); if prior_exist, if ~isempty(estim_params_.var_exo), @@ -341,6 +347,8 @@ if iload <=0, normaliz1 = [normaliz1; estim_params_.param_vals(:,7)]'; % normalize with prior standard deviation end % normaliz = max([normaliz; normaliz1]); + normaliz1(isinf(normaliz1)) = 1; + else normaliz1 = ones(1,nparam); end @@ -425,10 +433,24 @@ if iload <=0, if SampleSize > 1, - close(h) + close(h), + normTAU=std(TAU')'; + normLRE=std(LRE')'; + normGAM=std(GAM')'; + normaliz1=std(pdraws); + iter=0; + for ifile_index = 1:file_index, + load([IdentifDirectoryName '/' M_.fname '_identif_' int2str(ifile_index) '.mat'], 'stoH', 'stoJJ', 'stoLRE') + for irun=1:size(stoH,3), + iter=iter+1; + siJnorm(iter,:) = vnorm(stoJJ(:,:,irun)./repmat(normGAM,1,nparam)).*normaliz1; + siHnorm(iter,:) = vnorm(stoH(:,:,irun)./repmat(normTAU,1,nparam)).*normaliz1; + siLREnorm(iter,:) = vnorm(stoLRE(:,:,irun)./repmat(normLRE,1,nparam-offset)).*normaliz1(offset+1:end); + end + + end end - save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'pdraws', 'idemodel', 'idemoments', 'idelre', 'indJJ', 'indH', 'indLRE', ... 'siHmean', 'siJmean', 'siLREmean', 'derHmean', 'derJmean', 'derLREmean', 'TAU', 'GAM', 'LRE') else @@ -491,44 +513,44 @@ else close(h); end -if SampleSize>1, - siJmean = siJmean.*(ones(length(indJJ),1)*std(pdraws)); - siHmean = siHmean.*(ones(length(indH),1)*std(pdraws)); - siLREmean = siLREmean.*(ones(length(indLRE),1)*std(pdraws(:, offset+1:end ))); - - derJmean = derJmean.*(ones(length(indJJ),1)*std(pdraws)); - derHmean = derHmean.*(ones(length(indH),1)*std(pdraws)); - derLREmean = derLREmean.*(ones(length(indLRE),1)*std(pdraws(:, offset+1:end ))); - - derHmean = abs(derHmean./(max(siHmean')'*ones(1,size(pdraws,2)))); - derJmean = abs(derJmean./(max(siJmean')'*ones(1,size(pdraws,2)))); - derLREmean = abs(derLREmean./(max(siLREmean')'*ones(1,np))); - - siHmean = siHmean./(max(siHmean')'*ones(1,size(pdraws,2))); - siJmean = siJmean./(max(siJmean')'*ones(1,size(pdraws,2))); - siLREmean = siLREmean./(max(siLREmean')'*ones(1,np)); - - tstJmean = derJmean*0; - tstHmean = derHmean*0; - tstLREmean = derLREmean*0; - - if exist('OCTAVE_VERSION') - warning('off'), - else - warning off, - end - - for j=1:nparam, - indd = 1:length(siJmean(:,j)); - tstJmean(indd,j) = abs(derJmean(indd,j))./siJmean(indd,j); - indd = 1:length(siHmean(:,j)); - tstHmean(indd,j) = abs(derHmean(indd,j))./siHmean(indd,j); - if j>offset - indd = 1:length(siLREmean(:,j-offset)); - tstLREmean(indd,j-offset) = abs(derLREmean(indd,j-offset))./siLREmean(indd,j-offset); - end - end -end +% if SampleSize>1, +% siJmean = siJmean.*(ones(length(indJJ),1)*std(pdraws)); +% siHmean = siHmean.*(ones(length(indH),1)*std(pdraws)); +% siLREmean = siLREmean.*(ones(length(indLRE),1)*std(pdraws(:, offset+1:end ))); +% +% derJmean = derJmean.*(ones(length(indJJ),1)*std(pdraws)); +% derHmean = derHmean.*(ones(length(indH),1)*std(pdraws)); +% derLREmean = derLREmean.*(ones(length(indLRE),1)*std(pdraws(:, offset+1:end ))); +% +% derHmean = abs(derHmean./(max(siHmean')'*ones(1,size(pdraws,2)))); +% derJmean = abs(derJmean./(max(siJmean')'*ones(1,size(pdraws,2)))); +% derLREmean = abs(derLREmean./(max(siLREmean')'*ones(1,np))); +% +% siHmean = siHmean./(max(siHmean')'*ones(1,size(pdraws,2))); +% siJmean = siJmean./(max(siJmean')'*ones(1,size(pdraws,2))); +% siLREmean = siLREmean./(max(siLREmean')'*ones(1,np)); +% +% tstJmean = derJmean*0; +% tstHmean = derHmean*0; +% tstLREmean = derLREmean*0; +% +% if exist('OCTAVE_VERSION') +% warning('off'), +% else +% warning off, +% end +% +% for j=1:nparam, +% indd = 1:length(siJmean(:,j)); +% tstJmean(indd,j) = abs(derJmean(indd,j))./siJmean(indd,j); +% indd = 1:length(siHmean(:,j)); +% tstHmean(indd,j) = abs(derHmean(indd,j))./siHmean(indd,j); +% if j>offset +% indd = 1:length(siLREmean(:,j-offset)); +% tstLREmean(indd,j-offset) = abs(derLREmean(indd,j-offset))./siLREmean(indd,j-offset); +% end +% end +% end if nargout>3 && iload, @@ -547,141 +569,22 @@ end disp_identification(pdraws, idemodel, idemoments, name, advanced) - -if advanced, - figure('Name','Identification LRE model form'), - subplot(211) - if SampleSize > 1, - mmm = mean(siLREnorm); - else - mmm = (siLREnorm); - end - [ss, is] = sort(mmm); - if SampleSize ==1, - bar(siLREnorm(:,is)) - else - myboxplot(log(siLREnorm(:,is))) - end - % mmm = mean(siLREmean); - % [ss, is] = sort(mmm); - % myboxplot(siLREmean(:,is)) - % set(gca,'ylim',[0 1.05]) - set(gca,'xticklabel','') - for ip=1:np, - text(ip,-0.02,deblank(M_.param_names(indx(is(ip)),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') - end - title('Sensitivity in the LRE model') - - subplot(212) - if SampleSize>1, - mmm = mean(-idelre.Mco'); - else - mmm = (-idelre.Mco'); - end - [ss, is] = sort(mmm); - if SampleSize ==1, - bar(idelre.Mco(is,:)') - else - myboxplot(idelre.Mco(is,:)') - end - set(gca,'ylim',[0 1]) - set(gca,'xticklabel','') - for ip=1:np, - text(ip,-0.02,deblank(M_.param_names(indx(is(ip)),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') - end - title('Multicollinearity in the LRE model') - saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_LRE']) - eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_LRE']); - eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_LRE']); - if options_.nograph, close(gcf); end - - figure('Name','Identification in the model'), - % subplot(311) - % - % if SampleSize>1, - % mmm = mean(ide_strength_H); - % else - % mmm = (ide_strength_H); - % end - % [ss, is] = sort(mmm); - % if SampleSize>1, - % myboxplot(ide_strength_H(:,is)) - % else - % bar(ide_strength_H(:,is)) - % end - % % set(gca,'ylim',[0 1.05]) - % set(gca,'xticklabel','') - % dy = get(gca,'ylim'); - % % dy=dy(2)-dy(1); - % for ip=1:nparam, - % text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') - % end - % title('Identification strength in the model') - - subplot(211) - % mmm = mean(siHmean); - % [ss, is] = sort(mmm); - % myboxplot(siHmean(:,is)) - if SampleSize>1, - mmm = mean(siHnorm); - else - mmm = (siHnorm); - end - [ss, is] = sort(mmm); - if SampleSize>1, - myboxplot(log(siHnorm(:,is))) - else - bar(siHnorm(:,is)) - end - % set(gca,'ylim',[0 1.05]) - set(gca,'xticklabel','') - dy = get(gca,'ylim'); - % dy=dy(2)-dy(1); - for ip=1:nparam, - text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') - end - title('Sensitivity in the model') - - subplot(212) - if SampleSize>1, - mmm = mean(-idemodel.Mco'); - else - mmm = (-idemodel.Mco'); - end - % [ss, is] = sort(mmm); - if SampleSize>1, - myboxplot(idemodel.Mco(is,:)') - else - bar(idemodel.Mco(is,:)') - end - set(gca,'ylim',[0 1]) - set(gca,'xticklabel','') - for ip=1:nparam, - text(ip,-0.02,name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') - end - title('Multicollinearity in the model') - saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_model']) - eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_model']); - eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_model']); - if options_.nograph, close(gcf); end -end - figure('Name','Identification in the moments'), subplot(211) % mmm = mean(siHmean); % [ss, is] = sort(mmm); % myboxplot(siHmean(:,is)) -if SampleSize>1, - mmm = mean(ide_strength_J); -else +% if SampleSize>1, +% mmm = mean(ide_strength_J); +% else mmm = (ide_strength_J); -end +% end [ss, is] = sort(mmm); -if SampleSize>1, - myboxplot(log(ide_strength_J(:,is))) -else +% if SampleSize>1, +% myboxplot(log(ide_strength_J(:,is))) +% else bar(log([ide_strength_J(:,is)' ide_strength_J_prior(:,is)'])) -end +% end % set(gca,'ylim',[0 1.05]) set(gca,'xlim',[0 nparam+1]) set(gca,'xticklabel','') @@ -703,11 +606,11 @@ else mmm = (siJnorm); end % [ss, is] = sort(mmm); -if SampleSize > 1, - myboxplot(log(siJnorm(:,is))) -else - bar(siJnorm(:,is)) -end +% if SampleSize > 1, +% myboxplot(log(siJnorm(:,is))) +% else + bar(mmm(is)) +% end % set(gca,'ylim',[0 1.05]) set(gca,'xlim',[0 nparam+1]) set(gca,'xticklabel','') @@ -718,32 +621,62 @@ for ip=1:nparam, end title('Sensitivity in the moments') -% subplot(313) -% if SampleSize>1, -% mmm = mean(-idemoments.Mco'); -% else -% mmm = (-idemoments.Mco'); -% end -% % [ss, is] = sort(mmm); -% if SampleSize>1, -% myboxplot(idemoments.Mco(is,:)') -% else -% bar(idemoments.Mco(is,:)') -% end -% set(gca,'ylim',[0 1]) -% set(gca,'xticklabel','') -% for ip=1:nparam, -% text(ip,-0.02,name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') -% end -% title('Multicollinearity in the moments') -% saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_moments']) -% eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_moments']); -% eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_moments']); -% if options_.nograph, close(gcf); end +if advanced, + figure('Name','Identification in the model'), + subplot(211) + if SampleSize > 1, + mmm = mean(siLREnorm); + else + mmm = (siLREnorm); + end + mmm = [NaN(1, offset), mmm]; +% [ss, is] = sort(mmm); +% if SampleSize ==1, + bar(mmm(is)) +% else +% myboxplot(log(siLREnorm(:,is))) +% end + % mmm = mean(siLREmean); + % [ss, is] = sort(mmm); + % myboxplot(siLREmean(:,is)) + % set(gca,'ylim',[0 1.05]) + set(gca,'xticklabel','') +% for ip=1:np, +% text(ip,-0.02,deblank(M_.param_names(indx(is(ip)),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') + for ip=1:nparam, + text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') + end + title('Sensitivity in the LRE model') + + subplot(212) + % mmm = mean(siHmean); + % [ss, is] = sort(mmm); + % myboxplot(siHmean(:,is)) + if SampleSize>1, + mmm = mean(siHnorm); + else + mmm = (siHnorm); + end +% [ss, is] = sort(mmm); +% if SampleSize>1, +% myboxplot(log(siHnorm(:,is))) +% else + bar(mmm(is)) +% end + % set(gca,'ylim',[0 1.05]) + set(gca,'xticklabel','') + dy = get(gca,'ylim'); + % dy=dy(2)-dy(1); + for ip=1:nparam, + text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') + end + title('Sensitivity in the model') + + if options_.nograph, close(gcf); end -if SampleSize==1 && advanced, % identificaton patterns for j=1:size(cosnJ,2), + pax=NaN(nparam,nparam); fprintf('\n\n') disp(['Collinearity patterns with ', int2str(j) ,' parameter(s)']) fprintf('%-15s [%-*s] %10s\n','Parameter',(15+1)*j,' Expl. params ','cosn') @@ -753,7 +686,24 @@ if SampleSize==1 && advanced, namx=[namx ' ' sprintf('%-15s',name{pars{i,j}(in)})]; end fprintf('%-15s [%s] %10.3f\n',name{i},namx,cosnJ(i,j)) + pax(i,pars{i,j})=cosnJ(i,j); end + figure('name',['Collinearity patterns with ', int2str(j) ,' parameter(s)']), + imagesc(pax,[0 1]); + set(gca,'xticklabel','') + set(gca,'yticklabel','') + for ip=1:nparam, + 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 + colorbar; + ax=colormap; + ax(1,:)=[0.9 0.9 0.9]; + colormap(ax); + saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_collinearity_', int2str(j)]) + eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_collinearity_', int2str(j)]); + eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_collinearity_', int2str(j)]); + if options_.nograph, close(gcf); end end disp('') [U,S,V]=svd(siJ./normJ(:,ones(nparam,1)),0); @@ -779,8 +729,8 @@ if SampleSize==1 && advanced, % xlabel('SMALLEST singular values'), % end, else - bar(abs(V(:,j))), - Stit = S(j,j); + bar(abs(V(:,jj))), + Stit = S(jj,jj); % if j==8 || j==nparam, % xlabel('LARGEST singular values'), % end, @@ -803,10 +753,9 @@ if SampleSize==1 && advanced, eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_pattern_2']); eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_pattern_2']); end -end - -if advanced, - if SampleSize>1 + + if SampleSize>1, + options_.nograph=1; figure('Name','Condition Number'), subplot(221) hist(log10(idemodel.cond)) @@ -821,73 +770,26 @@ if advanced, eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_COND']); eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_COND']); if options_.nograph, close(gcf); end - end - - - pco = NaN(np,np); - for j=1:np, - if SampleSize>1 - pco(j+1:end,j) = mean(abs(squeeze(idelre.Pco(j+1:end,j,:))')); - else - pco(j+1:end,j) = abs(idelre.Pco(j+1:end,j))'; + ncut=floor(SampleSize/10*9); + [~,is]=sort(idelre.cond); + [proba, dproba] = stab_map_1(pdraws, is(1:ncut), is(ncut+1:end), 'HighestCondNumberLRE', 1, [], IdentifDirectoryName); + [~,is]=sort(idemodel.cond); + [proba, dproba] = stab_map_1(pdraws, is(1:ncut), is(ncut+1:end), 'HighestCondNumberModel', 1, [], IdentifDirectoryName); + [~,is]=sort(idemoments.cond); + [proba, dproba] = stab_map_1(pdraws, is(1:ncut), is(ncut+1:end), 'HighestCondNumberMoments', 1, [], IdentifDirectoryName); +% [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(pdraws, ibeh, inonbeh, ['HighestMultiCollinearity_',name{j}], 1, [], IdentifDirectoryName); +% end + [~,is]=sort(idemoments.Mco(j,:)); + [proba, dproba] = stab_map_1(pdraws, is(1:ncut), is(ncut+1:end), ['HighestMultiCollinearity_',name{j}], 1, [], IdentifDirectoryName); end end - figure('name','Pairwise correlations in the LRE model'), - imagesc(pco',[0 1]); - set(gca,'xticklabel','') - set(gca,'yticklabel','') - for ip=1:nparam, - 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 - colorbar; - saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_LRE']) - eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_LRE']); - eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_LRE']); - if options_.nograph, close(gcf); end - pco = NaN(nparam,nparam); - for j=1:nparam, - if SampleSize>1 - pco(j+1:end,j) = mean(abs(squeeze(idemodel.Pco(j+1:end,j,:))')); - else - pco(j+1:end,j) = abs(idemodel.Pco(j+1:end,j))'; - end - end - figure('name','Pairwise correlations in the model'), - imagesc(pco',[0 1]); - set(gca,'xticklabel','') - set(gca,'yticklabel','') - for ip=1:nparam, - 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 - colorbar; - saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_model']) - eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_model']); - eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_model']); - if options_.nograph, close(gcf); end - for j=1:nparam, - if SampleSize>1 - pco(j+1:end,j) = mean(abs(squeeze(idemoments.Pco(j+1:end,j,:))')); - else - pco(j+1:end,j) = abs(idemoments.Pco(j+1:end,j))'; - end - end - figure('name','Pairwise correlations in the moments'), - imagesc(pco',[0 1]); - set(gca,'xticklabel','') - set(gca,'yticklabel','') - for ip=1:nparam, - 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 - colorbar; - saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_moments']) - eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_moments']); - eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_moments']); - if options_.nograph, close(gcf); end end