plot_identification

 function plot_identification(params,idemoments,idehess,idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName)

 INPUTS
    o params             [array] parameter values for identification checks
    o idemoments         [structure] identification results for the moments
    o idehess            [structure] identification results for the Hessian
    o idemodel           [structure] identification results for the reduced form solution
    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 IdentifDirectoryName   [char] directory name
    
 OUTPUTS
    None
    
 SPECIAL REQUIREMENTS
    None

0001 function plot_identification(params,idemoments,idehess,idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName)
0002 % function plot_identification(params,idemoments,idehess,idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName)
0003 %
0004 % INPUTS
0005 %    o params             [array] parameter values for identification checks
0006 %    o idemoments         [structure] identification results for the moments
0007 %    o idehess            [structure] identification results for the Hessian
0008 %    o idemodel           [structure] identification results for the reduced form solution
0009 %    o idelre             [structure] identification results for the LRE model
0010 %    o advanced           [integer] flag for advanced identification checks
0011 %    o tittxt             [char] name of the results to plot
0012 %    o name               [char] list of names
0013 %    o IdentifDirectoryName   [char] directory name
0014 %
0015 % OUTPUTS
0016 %    None
0017 %
0018 % SPECIAL REQUIREMENTS
0019 %    None
0020 
0021 % Copyright (C) 2008-2011 Dynare Team
0022 %
0023 % This file is part of Dynare.
0024 %
0025 % Dynare is free software: you can redistribute it and/or modify
0026 % it under the terms of the GNU General Public License as published by
0027 % the Free Software Foundation, either version 3 of the License, or
0028 % (at your option) any later version.
0029 %
0030 % Dynare is distributed in the hope that it will be useful,
0031 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0032 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0033 % GNU General Public License for more details.
0034 %
0035 % You should have received a copy of the GNU General Public License
0036 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0037 
0038 global M_ options_
0039 
0040 
0041 [SampleSize, nparam]=size(params);
0042 siJnorm = idemoments.siJnorm;
0043 siHnorm = idemodel.siHnorm;
0044 siLREnorm = idelre.siLREnorm;
0045 
0046 % if prior_exist,
0047 %     tittxt = 'Prior mean - ';
0048 % else
0049 %     tittxt = '';
0050 % end
0051 tittxt1=regexprep(tittxt, ' ', '_');
0052 tittxt1=strrep(tittxt1, '.', '');
0053 if SampleSize == 1,
0054     siJ = idemoments.siJ;
0055     hh = dyn_figure(options_,'Name',[tittxt, ' - Identification using info from observables']);
0056     subplot(211)
0057     mmm = (idehess.ide_strength_J);
0058     [ss, is] = sort(mmm);
0059     bar(log([idehess.ide_strength_J(:,is)' idehess.ide_strength_J_prior(:,is)']))
0060     set(gca,'xlim',[0 nparam+1])
0061     set(gca,'xticklabel','')
0062     dy = get(gca,'ylim');
0063     for ip=1:nparam,
0064         text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
0065     end
0066     legend('relative to param value','relative to prior std','Location','Best')
0067     if  idehess.flag_score,
0068         title('Identification strength with asymptotic Information matrix (log-scale)')
0069     else
0070         title('Identification strength with moments Information matrix (log-scale)')
0071     end
0072     
0073     subplot(212)
0074     bar(log([idehess.deltaM(is) idehess.deltaM_prior(is)]))
0075     set(gca,'xlim',[0 nparam+1])
0076     set(gca,'xticklabel','')
0077     dy = get(gca,'ylim');
0078     for ip=1:nparam,
0079         text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
0080     end
0081     legend('relative to param value','relative to prior std','Location','Best')
0082     if  idehess.flag_score,
0083         title('Sensitivity component with asymptotic Information matrix (log-scale)')
0084     else
0085         title('Sensitivity component with moments Information matrix (log-scale)')
0086     end
0087     dyn_saveas(hh,[IdentifDirectoryName '/' M_.fname '_ident_strength_' tittxt1],options_);
0088     
0089     if advanced,
0090         disp(' ')
0091         disp('Press ENTER to display advanced diagnostics'), pause(5),
0092         hh = dyn_figure(options_,'Name',[tittxt, ' - Sensitivity plot']);
0093         subplot(211)
0094         mmm = (siJnorm)'./max(siJnorm);
0095         mmm1 = (siHnorm)'./max(siHnorm);
0096         mmm=[mmm mmm1];
0097         mmm1 = (siLREnorm)'./max(siLREnorm);
0098         offset=length(siHnorm)-length(siLREnorm);
0099         mmm1 = [NaN(offset,1); mmm1];
0100         mmm=[mmm mmm1];
0101         
0102         bar(log(mmm(is,:).*100))
0103         set(gca,'xlim',[0 nparam+1])
0104         set(gca,'xticklabel','')
0105         dy = get(gca,'ylim');
0106         for ip=1:nparam,
0107             text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
0108         end
0109         title('Sensitivity bars using derivatives (log-scale)')
0110         dyn_saveas(hh,[IdentifDirectoryName '/' M_.fname '_sensitivity_' tittxt1 ],options_);
0111         
0112         % identificaton patterns
0113         for  j=1:size(idemoments.cosnJ,2),
0114             pax=NaN(nparam,nparam);
0115             fprintf('\n')
0116             disp(['Collinearity patterns with ', int2str(j) ,' parameter(s)'])
0117             fprintf('%-15s [%-*s] %10s\n','Parameter',(15+1)*j,' Expl. params ','cosn')
0118             for i=1:nparam,
0119                 namx='';
0120                 for in=1:j,
0121                     dumpindx = idemoments.pars{i,j}(in);
0122                     if isnan(dumpindx),
0123                         namx=[namx ' ' sprintf('%-15s','--')];
0124                     else
0125                         namx=[namx ' ' sprintf('%-15s',name{dumpindx})];
0126                         pax(i,dumpindx)=idemoments.cosnJ(i,j);
0127                     end
0128                 end
0129                 fprintf('%-15s [%s] %10.3f\n',name{i},namx,idemoments.cosnJ(i,j))
0130             end
0131             hh = dyn_figure(options_,'Name',[tittxt,' - Collinearity patterns with ', int2str(j) ,' parameter(s)']);
0132             imagesc(pax,[0 1]);
0133             set(gca,'xticklabel','')
0134             set(gca,'yticklabel','')
0135             for ip=1:nparam,
0136                 text(ip,(0.5),name{ip},'rotation',90,'HorizontalAlignment','left','interpreter','none')
0137                 text(0.5,ip,name{ip},'rotation',0,'HorizontalAlignment','right','interpreter','none')
0138             end
0139             colorbar;
0140             ax=colormap;
0141             ax(1,:)=[0.9 0.9 0.9];
0142             colormap(ax);
0143             if nparam>10,
0144                 set(gca,'xtick',(5:5:nparam))
0145                 set(gca,'ytick',(5:5:nparam))
0146             end
0147             set(gca,'xgrid','on')
0148             set(gca,'ygrid','on')
0149             xlabel([tittxt,' - Collinearity patterns with ', int2str(j) ,' parameter(s)'],'interpreter','none')
0150             dyn_saveas(hh,[ IdentifDirectoryName '/' M_.fname '_ident_collinearity_' tittxt1 '_' int2str(j) ],options_);
0151         end
0152         disp('')
0153         if idehess.flag_score,
0154             [U,S,V]=svd(idehess.AHess,0);
0155             S=diag(S);
0156             if nparam<5,
0157                 f1 = dyn_figure(options_,'Name',[tittxt,' - Identification patterns (Information matrix)']);
0158             else
0159                 f1 = dyn_figure(options_,'Name',[tittxt,' - Identification patterns (Information matrix): SMALLEST SV']);
0160                 f2 = dyn_figure(options_,'Name',[tittxt,' - Identification patterns (Information matrix): HIGHEST SV']);
0161             end
0162         else
0163             S = idemoments.S;
0164             V = idemoments.V;
0165             if nparam<5,
0166                 f1 = dyn_figure(options_,'Name',[tittxt,' - Identification patterns (moments)']);
0167             else
0168                 f1 = dyn_figure(options_,'Name',[tittxt,' - Identification patterns (moments): SMALLEST SV']);
0169                 f2 = dyn_figure(options_,'Name',[tittxt,' - Identification patterns (moments): HIGHEST SV']);
0170             end
0171         end
0172         for j=1:min(nparam,8),
0173             if j<5,
0174                 figure(f1),
0175                 jj=j;
0176             else
0177                 figure(f2),
0178                 jj=j-4;
0179             end
0180             subplot(4,1,jj),
0181             if j<5
0182                 bar(abs(V(:,end-j+1))),
0183                 Stit = S(end-j+1);
0184             else
0185                 bar(abs(V(:,jj))),
0186                 Stit = S(jj);
0187             end
0188             set(gca,'xticklabel','')
0189             if j==4 || j==nparam || j==8,
0190                 for ip=1:nparam,
0191                     text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
0192                 end
0193             end
0194             title(['Singular value ',num2str(Stit)])
0195         end
0196         dyn_saveas(f1,[  IdentifDirectoryName '/' M_.fname '_ident_pattern_' tittxt1 '_1' ],options_);
0197         if nparam>4,
0198             dyn_saveas(f2,[  IdentifDirectoryName '/' M_.fname '_ident_pattern_' tittxt1 '_2' ],options_);
0199         end
0200     end
0201     
0202 else
0203     hh = dyn_figure(options_,'Name',['MC sensitivities']);
0204     subplot(211)
0205     mmm = (idehess.ide_strength_J);
0206     [ss, is] = sort(mmm);
0207     mmm = mean(siJnorm)';
0208     mmm = mmm./max(mmm);
0209     if advanced,
0210         mmm1 = mean(siHnorm)';
0211         mmm=[mmm mmm1./max(mmm1)];
0212         mmm1 = mean(siLREnorm)';
0213         offset=size(siHnorm,2)-size(siLREnorm,2);
0214         mmm1 = [NaN(offset,1); mmm1./max(mmm1)];
0215         mmm=[mmm mmm1];
0216     end        
0217         
0218     bar(mmm(is,:))
0219     set(gca,'xlim',[0 nparam+1])
0220     set(gca,'xticklabel','')
0221     dy = get(gca,'ylim');
0222     for ip=1:nparam,
0223         text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
0224     end
0225     if advanced,
0226         legend('Moments','Model','LRE model','Location','Best')
0227     end
0228     title('MC mean of sensitivity measures')
0229     dyn_saveas(hh,[ IdentifDirectoryName '/' M_.fname '_MC_sensitivity' ],options_);
0230     if advanced,
0231         disp(' ')
0232         disp('Press ENTER to display advanced diagnostics'), pause(5),
0233 %         options_.nograph=1;
0234         hh = dyn_figure(options_,'Name','MC Condition Number');
0235         subplot(221)
0236         hist(log10(idemodel.cond))
0237         title('log10 of Condition number in the model')
0238         subplot(222)
0239         hist(log10(idemoments.cond))
0240         title('log10 of Condition number in the moments')
0241         subplot(223)
0242         hist(log10(idelre.cond))
0243         title('log10 of Condition number in the LRE model')
0244         dyn_saveas(hh,[IdentifDirectoryName '/' M_.fname '_ident_COND' ],options_);
0245         ncut=floor(SampleSize/10*9);
0246         [dum,is]=sort(idelre.cond);
0247         [proba, dproba] = stab_map_1(params, is(1:ncut), is(ncut+1:end), 'MC_HighestCondNumberLRE', 1, [], IdentifDirectoryName, 0.1);
0248         [dum,is]=sort(idemodel.cond);
0249         [proba, dproba] = stab_map_1(params, is(1:ncut), is(ncut+1:end), 'MC_HighestCondNumberModel', 1, [], IdentifDirectoryName, 0.1);
0250         [dum,is]=sort(idemoments.cond);
0251         [proba, dproba] = stab_map_1(params, is(1:ncut), is(ncut+1:end), 'MC_HighestCondNumberMoments', 1, [], IdentifDirectoryName, 0.1);
0252 %         [proba, dproba] = stab_map_1(idemoments.Mco', is(1:ncut), is(ncut+1:end), 'HighestCondNumberMoments_vs_Mco', 1, [], IdentifDirectoryName);
0253 %         for j=1:nparam,
0254 % %             ibeh=find(idemoments.Mco(j,:)<0.9);
0255 % %             inonbeh=find(idemoments.Mco(j,:)>=0.9);
0256 % %             if ~isempty(ibeh) && ~isempty(inonbeh)
0257 % %                 [proba, dproba] = stab_map_1(params, ibeh, inonbeh, ['HighestMultiCollinearity_',name{j}], 1, [], IdentifDirectoryName);
0258 % %             end
0259 %             [~,is]=sort(idemoments.Mco(:,j));
0260 %             [proba, dproba] = stab_map_1(params, is(1:ncut), is(ncut+1:end), ['MC_HighestMultiCollinearity_',name{j}], 1, [], IdentifDirectoryName, 0.15);
0261 %         end
0262 
0263         if nparam<5,
0264             f1 = dyn_figure(options_,'Name',[tittxt,' - MC Identification patterns (moments): HIGHEST SV'])
0265         else
0266             f1 = dyn_figure(options_,'Name',[tittxt,' - MC Identification patterns (moments): SMALLEST SV']);
0267             f2 = dyn_figure(options_,'Name',[tittxt,' - MC Identification patterns (moments): HIGHEST SV']);
0268         end
0269         nplots=min(nparam,8);
0270         if nplots>4,
0271             nsubplo=ceil(nplots/2);
0272         else
0273             nsubplo=nplots;
0274         end
0275         for j=1:nplots,
0276             if (nparam>4 && j<=ceil(nplots/2)) || nparam<5,
0277                 figure(f1),
0278                 jj=j;
0279                 VVV=squeeze(abs(idemoments.V(:,:,end-j+1)));
0280                 SSS = idemoments.S(:,end-j+1);
0281             else
0282                 figure(f2),
0283                 jj=j-ceil(nplots/2);
0284                 VVV=squeeze(abs(idemoments.V(:,:,jj)));
0285                 SSS = idemoments.S(:,jj);
0286             end
0287             subplot(nsubplo,1,jj),
0288             for i=1:nparam,
0289                 [post_mean, post_median(:,i), post_var, hpd_interval(i,:), post_deciles] = posterior_moments(VVV(:,i),0,0.9);
0290             end
0291             bar(post_median)
0292             hold on, plot(hpd_interval,'--*r'),
0293             Stit=mean(SSS);
0294 
0295             set(gca,'xticklabel','')
0296             if j==4 || j==nparam || j==8,
0297                 for ip=1:nparam,
0298                     text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
0299                 end
0300             end
0301             title(['MEAN Singular value ',num2str(Stit)])
0302         end
0303         dyn_saveas(f1,[IdentifDirectoryName '/' M_.fname '_MC_ident_pattern_1' ],options_);
0304         if nparam>4,
0305             dyn_saveas(f2,[  IdentifDirectoryName '/' M_.fname '_MC_ident_pattern_2' ],options_);
0306         end
0307     end
0308 end
0309 
0310 % disp_identification(params, idemodel, idemoments, name)