diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m
index 55c1186b1..5ef97a543 100644
--- a/matlab/dynare_identification.m
+++ b/matlab/dynare_identification.m
@@ -30,14 +30,27 @@ end
 options_ident = set_default_option(options_ident,'load_ident_files',0);
 options_ident = set_default_option(options_ident,'useautocorr',0);
 options_ident = set_default_option(options_ident,'ar',3);
-options_ident = set_default_option(options_ident,'prior_mc',2000);
+options_ident = set_default_option(options_ident,'prior_mc',1);
 options_ident = set_default_option(options_ident,'prior_range',0);
+options_ident = set_default_option(options_ident,'periods',300);
+options_ident = set_default_option(options_ident,'replic',100);
+options_ident = set_default_option(options_ident,'advanced',0);
 if nargin==2,
     options_ident.prior_mc=size(pdraws0,1);
 end
+if isempty(estim_params_),
+    options_ident.prior_mc=1;
+    options_ident.prior_range=0;
+    prior_exist=0;
+else
+    prior_exist=1;
+end    
 
 iload = options_ident.load_ident_files;
+advanced = options_ident.advanced;
 nlags = options_ident.ar;
+periods = options_ident.periods;
+replic = options_ident.replic;
 useautocorr = options_ident.useautocorr;
 options_.ar=nlags;
 options_.prior_mc = options_ident.prior_mc;
@@ -68,33 +81,60 @@ if ~(exist('sylvester3mr','file')==2),
 end
 
 IdentifDirectoryName = CheckPath('identification');
+if prior_exist,
 
-indx = estim_params_.param_vals(:,1);
-indexo=[];
-if ~isempty(estim_params_.var_exo)
-    indexo = estim_params_.var_exo(:,1);
+    indx = estim_params_.param_vals(:,1);
+    indexo=[];
+    if ~isempty(estim_params_.var_exo)
+        indexo = estim_params_.var_exo(:,1);
+    end
+
+    nparam = length(bayestopt_.name);
+    np = estim_params_.np;
+    name = bayestopt_.name;
+
+    offset = estim_params_.nvx;
+    offset = offset + estim_params_.nvn;
+    offset = offset + estim_params_.ncx;
+    offset = offset + estim_params_.ncn;
+else
+    indx = [1:M_.param_nbr];
+    indexo = [1:M_.exo_nbr];
+    offset = M_.exo_nbr;
+    np = M_.param_nbr;
+    nparam = np+offset;
+    name = [cellstr(M_.exo_names); cellstr(M_.param_names)];
 end
 
-nparam = length(bayestopt_.name);
-
 MaxNumberOfBytes=options_.MaxNumberOfBytes;
 
+disp(' ')
+disp(['==== Identification analysis ====' ]),
+disp(' ')
 
 if iload <=0, 
     
     iteration = 0;
     burnin_iteration = 0;
+    if SampleSize==1,
+        BurninSampleSize=0;
+    else
+        BurninSampleSize=50;
+    end
     loop_indx = 0;
     file_index = 0;
     run_index = 0;
 
+    if SampleSize > 1,
     h = waitbar(0,'Monte Carlo identification checks ...');
+    end
     [I,J]=find(M_.lead_lag_incidence');
 
     while iteration < SampleSize,
         loop_indx = loop_indx+1;
+        if prior_exist,
         if nargin==2,
-            if burnin_iteration>=50,
+            if burnin_iteration>=BurninSampleSize,
                 params = pdraws0(iteration+1,:);
             else
                 params = pdraws0(burnin_iteration+1,:);
@@ -103,6 +143,9 @@ if iload <=0,
             params = prior_draw();
         end
         set_all_parameters(params);
+        else
+            params = [sqrt(diag(M_.Sigma_e))', M_.params'];
+        end
         [A,B,ys,info]=dynare_resolve;
 
         
@@ -117,7 +160,7 @@ if iload <=0,
             yy0=oo_.dr.ys(I);    
             [residual, g1 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', M_.params,1);    
 
-            if burnin_iteration<50,
+            if burnin_iteration<BurninSampleSize,
                 burnin_iteration = burnin_iteration + 1;
                 pdraws(burnin_iteration,:) = params;
                 TAU(:,burnin_iteration)=tau;
@@ -143,11 +186,11 @@ if iload <=0,
                     dum = [dum; vec(gam{j+1})];
                 end
                 GAM(:,burnin_iteration)=[oo_.dr.ys(bayestopt_.mfys); dum];
-                warning warning_old_state;
+%                 warning warning_old_state;
             else
                 iteration = iteration + 1;
                 run_index = run_index + 1;
-                if iteration==1,
+                if iteration==1 & BurninSampleSize,
                     indJJ = (find(std(GAM')>1.e-8));
                     indH = (find(std(TAU')>1.e-8));
                     indLRE = (find(std(LRE')>1.e-8));
@@ -163,10 +206,23 @@ if iload <=0,
             end
 
             if iteration,
+                [JJ, H, gam, gp] = getJJ(A, B, M_,oo0,options_,0,indx,indexo,bayestopt_.mf2,nlags,useautocorr);      
+                if BurninSampleSize == 0,
+                    indJJ = (find(max(abs(JJ'))>1.e-8));
+                    indH = (find(max(abs(H'))>1.e-8));
+                    indLRE = (find(max(abs(gp'))>1.e-8));
+                    TAU = zeros(length(indH),SampleSize);
+                    GAM = zeros(length(indJJ),SampleSize);
+                    LRE = zeros(length(indLRE),SampleSize);
+                    MAX_tau   = min(SampleSize,ceil(MaxNumberOfBytes/(length(indH)*nparam)/8));
+                    MAX_gam   = min(SampleSize,ceil(MaxNumberOfBytes/(length(indJJ)*nparam)/8));
+                    stoH = zeros([length(indH),nparam,MAX_tau]);
+                    stoJJ = zeros([length(indJJ),nparam,MAX_tau]);
+                    delete([IdentifDirectoryName '/' M_.fname '_identif_*.mat'])
+                end
                 TAU(:,iteration)=tau(indH);
                 vg1 = [oo_.dr.ys(oo_.dr.order_var); vec(g1)];
                 LRE(:,iteration)=vg1(indLRE);
-                [JJ, H, gam, gp] = getJJ(A, B, M_,oo0,options_,0,indx,indexo,bayestopt_.mf2,nlags,useautocorr);      
                 GAM(:,iteration)=gam(indJJ);
                 stoLRE(:,:,run_index) = gp(indLRE,:);
                 stoH(:,:,run_index) = H(indH,:);
@@ -177,12 +233,16 @@ if iload <=0,
                 % use prior uncertainty
                 siJ = (JJ(indJJ,:));
                 siH = (H(indH,:));
+
                 siLRE = (gp(indLRE,:));
                 %       siJ = abs(JJ(indJJ,:).*(ones(length(indJJ),1)*bayestopt_.p2'));
                 %       siH = abs(H(indH,:).*(ones(length(indH),1)*bayestopt_.p2'));
                 %       siJ = abs(JJ(indJJ,:).*(1./mGAM'*bayestopt_.p2'));
                 %       siH = abs(H(indH,:).*(1./mTAU'*bayestopt_.p2'));
 
+                siJnorm(iteration,:) = vnorm(siJ./repmat(GAM(:,iteration),1,nparam)).*params;        
+                siHnorm(iteration,:) = vnorm(siH./repmat(TAU(:,iteration),1,nparam)).*params;        
+                siLREnorm(iteration,:) = vnorm(siLRE./repmat(LRE(:,iteration),1,nparam-offset)).*params(offset+1:end);        
                 if iteration ==1,
                     siJmean = abs(siJ)./SampleSize;        
                     siHmean = abs(siH)./SampleSize;        
@@ -210,8 +270,36 @@ if iload <=0,
                  idemodel.cond(iteration), idemoments.cond(iteration), idelre.cond(iteration), ...        
                  idemodel.ee(:,:,iteration), idemoments.ee(:,:,iteration), idelre.ee(:,:,iteration), ...        
                  idemodel.ind(:,iteration), idemoments.ind(:,iteration), ...
-                 idemodel.indno{iteration}, idemoments.indno{iteration}] = ...
+                idemodel.indno{iteration}, idemoments.indno{iteration}, ...
+                idemodel.ino(iteration), idemoments.ino(iteration)] = ...
                  identification_checks(H(indH,:)./normH(:,ones(nparam,1)),JJ(indJJ,:)./normJ(:,ones(nparam,1)), gp(indLRE,:)./normLRE(:,ones(size(gp,2),1)), bayestopt_);      
+%                  identification_checks(H(indH,:),JJ(indJJ,:), gp(indLRE,:), bayestopt_);      
+                indok = find(max(idemoments.indno{iteration},[],1)==0);
+                ide_strength_J(iteration,:)=NaN(1,nparam);
+                if iteration ==1,
+                    cmm = simulated_moment_uncertainty(indJJ, periods, replic);
+                end,
+%                 Jinv=(siJ(:,indok)'*siJ(:,indok))\siJ(:,indok)';
+%                 MIM=inv(Jinv*cmm*Jinv');
+                MIM=siJ(:,indok)'*inv(cmm)*siJ(:,indok);
+                deltaM = sqrt(diag(MIM));
+                tildaM = MIM./((deltaM)*(deltaM'));
+                rhoM=sqrt(1-1./diag(inv(tildaM)));
+                deltaM = deltaM.*params(indok)';
+                ide_strength_J(iteration,indok) = (1./[sqrt(diag(inv(MIM)))./params(indok)']);
+
+%                 indok = find(max(idemodel.indno{iteration},[],1)==0);
+%                 ide_strength_H(iteration,:)=zeros(1,nparam);
+%                 mim=inv(siH(:,indok)'*siH(:,indok))*siH(:,indok)';
+% %                 mim=mim*diag(GAM(:,iteration))*mim';
+% %                 MIM=inv(mim);
+%                 mim=mim.*repmat(TAU(:,iteration),1,length(indok))';
+%                 MIM=inv(mim*mim');
+%                 deltaM = sqrt(diag(MIM));
+%                 tildaM = MIM./((deltaM)*(deltaM'));
+%                 rhoM=sqrt(1-1./diag(inv(tildaM)));
+%                 deltaM = deltaM.*params(indok)';
+%                 ide_strength_H(iteration,indok) = (1./[sqrt(diag(inv(MIM)))./params(indok)']);
                 if run_index==MAX_tau | iteration==SampleSize,
                     file_index = file_index + 1;
                     if run_index<MAX_tau,
@@ -224,13 +312,17 @@ if iload <=0,
                     
                 end
 
+                if SampleSize > 1,
                 waitbar(iteration/SampleSize,h,['MC Identification checks ',int2str(iteration),'/',int2str(SampleSize)])
+                end
             end
         end
     end
 
 
+    if SampleSize > 1,
     close(h)
+    end
 
 
     save([IdentifDirectoryName '/' M_.fname '_identif'], 'pdraws', 'idemodel', 'idemoments', 'idelre', 'indJJ', 'indH', 'indLRE', ...  
@@ -263,7 +355,8 @@ else
                 idemodel.cond(iteration), idemoments.cond(iteration), idelre.cond(iteration), ...        
                 idemodel.ee(:,:,iteration), idemoments.ee(:,:,iteration), idelre.ee(:,:,iteration), ...        
                 idemodel.ind(:,iteration), idemoments.ind(:,iteration), ...
-                idemodel.indno{iteration}, idemoments.indno{iteration}] = ...
+                idemodel.indno{iteration}, idemoments.indno{iteration}, ...
+                idemodel.ino(iteration), idemoments.ino(iteration)] = ...
                 identification_checks(stoH(:,:,index)./normH(:,ones(nparam,1)), ...
                 stoJJ(:,:,index)./normJ(:,ones(nparam,1)), ...
                 stoLRE(:,:,index)./normLRE(:,ones(size(stoLRE,2),1)), bayestopt_);      
@@ -273,13 +366,28 @@ else
     close(h);
     save([IdentifDirectoryName '/' M_.fname '_identif'], 'idemodel', 'idemoments', 'idelre', '-append')
     end
+    iteration = 0;
+    h = waitbar(0,'Monte Carlo identification checks ...');
+    for file_index=1:length(identFiles)
+        load([IdentifDirectoryName '/' M_.fname '_identif_' int2str(file_index)], 'stoH', 'stoJJ', 'stoLRE')      
+        for index=1:size(stoH,3),
+            iteration = iteration+1;
+            fobj(iteration)=sum(((GAM(:,iteration)-GAM(:,1))).^2);
+            fobjH(iteration)=sum(((TAU(:,iteration)-TAU(:,1))).^2);
+            fobjR(iteration)=sum(((GAM(:,iteration)./GAM(:,1)-1)).^2);
+            fobjHR(iteration)=sum(((TAU(:,iteration)./TAU(:,1)-1)).^2);
+            FOC(iteration,:) = (GAM(:,iteration)-GAM(:,1))'*stoJJ(:,:,index);
+            FOCH(iteration,:) = (TAU(:,iteration)-TAU(:,1))'*stoH(:,:,index);
+            FOCR(iteration,:) = ((GAM(:,iteration)./GAM(:,1)-1)./GAM(:,1))'*stoJJ(:,:,index);
+            FOCHR(iteration,:) = ((TAU(:,iteration)./TAU(:,1)-1)./TAU(:,1))'*stoH(:,:,index);
+
+            waitbar(iteration/SampleSize,h,['MC Identification checks ',int2str(iteration),'/',int2str(SampleSize)])
+        end
+    end
+    close(h);
 end  
 
-offset = estim_params_.nvx;
-offset = offset + estim_params_.nvn;
-offset = offset + estim_params_.ncx;
-offset = offset + estim_params_.ncn;
-
+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 )));
@@ -290,15 +398,16 @@ 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,estim_params_.np)));
+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,estim_params_.np));
+siLREmean = siLREmean./(max(siLREmean')'*ones(1,np));
 
 tstJmean = derJmean*0;
 tstHmean = derHmean*0;
 tstLREmean = derLREmean*0;
+
 if exist('OCTAVE_VERSION')
     warning('off', 'Octave:divide-by-zero')
 else
@@ -315,6 +424,7 @@ if j>offset
 tstLREmean(indd,j-offset) = abs(derLREmean(indd,j-offset))./siLREmean(indd,j-offset);
 end  
 end
+end
 
 
 if nargout>3 & iload,
@@ -382,7 +492,7 @@ end
 % siHPCA2 = siHPCA2./max(siHPCA2);
 
 
-disp_identification(pdraws, idemodel, idemoments)
+disp_identification(pdraws, idemodel, idemoments, name, advanced)
 
 % figure,
 % % myboxplot(siPCA(1:(max(find(cumsum(latent)./length(indJJ)<0.99))+1),:))
@@ -392,7 +502,7 @@ disp_identification(pdraws, idemodel, idemoments)
 % set(gca,'xticklabel','')
 % set(gca,'xlim',[0.5 nparam+0.5])
 % for ip=1:nparam,
-%   text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+%   text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
 % end
 % title('Sensitivity in TAU''s PCA')
 % 
@@ -403,7 +513,7 @@ disp_identification(pdraws, idemodel, idemoments)
 % set(gca,'xticklabel','')
 % set(gca,'xlim',[0.5 nparam+0.5])
 % for ip=1:nparam,
-%   text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+%   text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
 % end
 % title('Sensitivity in standardized TAU''s PCA')
 % 
@@ -415,7 +525,7 @@ disp_identification(pdraws, idemodel, idemoments)
 % set(gca,'xticklabel','')
 % set(gca,'xlim',[0.5 nparam+0.5])
 % for ip=1:nparam,
-%   text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+%   text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
 % end
 % title('Sensitivity in moments'' PCA')
 % 
@@ -426,30 +536,50 @@ disp_identification(pdraws, idemodel, idemoments)
 % set(gca,'xticklabel','')
 % set(gca,'xlim',[0.5 nparam+0.5])
 % for ip=1:nparam,
-%   text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+%   text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
 % end
 % title('Sensitivity in standardized moments'' PCA')
 
+if advanced,
 figure('Name','Identification LRE model form'),
 subplot(211)
-mmm = mean(siLREmean);
+if SampleSize > 1,
+mmm = mean(siLREnorm);
+else
+mmm = (siLREnorm);
+end
 [ss, is] = sort(mmm);
-myboxplot(siLREmean(:,is))
-set(gca,'ylim',[0 1.05])
+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:estim_params_.np,
-  text(ip,-0.02,deblank(M_.param_names(estim_params_.param_vals(is(ip),1),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
+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);
-myboxplot(idelre.Mco(is,:)')
+if SampleSize ==1,
+    bar(idelre.Mco(is,:)')
+else
+    myboxplot(idelre.Mco(is,:)')
+end
 set(gca,'ylim',[0 1])
 set(gca,'xticklabel','')
-for ip=1:estim_params_.np,
-  text(ip,-0.02,deblank(M_.param_names(estim_params_.param_vals(is(ip),1),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
+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'])
@@ -458,60 +588,176 @@ 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);
+% mmm = mean(siHmean);
+% [ss, is] = sort(mmm);
+% myboxplot(siHmean(:,is))
+if SampleSize>1,
+mmm = mean(siHnorm);
+else
+mmm = (siHnorm);
+end
 [ss, is] = sort(mmm);
-myboxplot(siHmean(:,is))
-set(gca,'ylim',[0 1.05])
+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,-0.02,bayestopt_.name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+    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');
-[ss, is] = sort(mmm);
+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,bayestopt_.name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+    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(siJmean);
+% mmm = mean(siHmean);
+% [ss, is] = sort(mmm);
+% myboxplot(siHmean(:,is))
+if SampleSize>1,
+mmm = mean(ide_strength_J);
+else
+mmm = (ide_strength_J);
+end
 [ss, is] = sort(mmm);
-myboxplot(siJmean(:,is))
-set(gca,'ylim',[0 1.05])
+if SampleSize>1,
+myboxplot(log(ide_strength_J(:,is)))
+else
+bar(ide_strength_J(:,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,-0.02,bayestopt_.name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+    text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+end
+title('Identification strength in the moments')
+
+subplot(212)
+% mmm = mean(siJmean);
+% [ss, is] = sort(mmm);
+% myboxplot(siJmean(:,is))
+if SampleSize > 1,
+mmm = mean(siJnorm);
+else
+mmm = (siJnorm);
+end
+% [ss, is] = sort(mmm);
+if SampleSize > 1,
+myboxplot(log(siJnorm(:,is)))
+else
+bar(siJnorm(:,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 moments')
 
-subplot(212)
-mmm = mean(-idemoments.Mco');
-[ss, is] = sort(mmm);
-myboxplot(idemoments.Mco(is,:)')
-set(gca,'ylim',[0 1])
-set(gca,'xticklabel','')
-for ip=1:nparam,
-    text(ip,-0.02,bayestopt_.name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+% 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 SampleSize==1,
+    % identificaton patterns
+    [U,S,V]=svd(siJ./normJ(:,ones(nparam,1)),0);
+    figure('name','Identification patterns (moments)'),
+    for j=1:min(nparam,8),
+        subplot(2,4,j),
+        if j<5
+            bar(abs(V(:,end-j+1))),
+            Stit = S(end-j+1,end-j+1);
+            if j==1, ylabel('SMALLEST singular values'), end,
+        else
+            bar(abs(V(:,j))),
+            Stit = S(j,j);
+            if j==5, ylabel('LARGEST singular values'), end,
+        end
+        set(gca,'xticklabel','')
+        for ip=1:nparam,
+            text(ip,-0.02,name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+        end
+        title(['Singular value ',num2str(Stit)])
+    end
+    saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_pattern'])
+    eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_pattern']);
+    eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_pattern']);
 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,
+if SampleSize>1
 figure('Name','Condition Number'),
 subplot(221)
 hist(log10(idemodel.cond))
@@ -526,98 +772,178 @@ saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_COND'])
 eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_COND']);
 eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_COND']);
 if options_.nograph, close(gcf); end
-ifig=0;
-nbox = min(estim_params_.np-1,12);
-for j=1:estim_params_.np,  
-    if mod(j,12)==1,    
-        ifig = ifig+1;    
-        figure('name','Partial correlations in the LRE model'),    
-        iplo=0;  
-    end
-    iplo=iplo+1;    
-    mmm = mean(squeeze(idelre.Pco(:,j,:))');    
-    [sss, immm] = sort(-mmm);    
-    subplot(3,4,iplo), 
-    if nbox==1,
-        myboxplot(squeeze(idelre.Pco(immm(2:nbox+1),j,:))),
-    else
-        myboxplot(squeeze(idelre.Pco(immm(2:nbox+1),j,:))'),
-    end
-    set(gca,'ylim',[0 1])
-    set(gca,'xticklabel','')
-    for ip=1:nbox, %estim_params_.np,  
-        text(ip,-0.02,deblank(M_.param_names(estim_params_.param_vals(immm(ip+1),1),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
-    end
-    title(deblank(M_.param_names(estim_params_.param_vals(j,1),:))),
-    if j==estim_params_.np | mod(j,12)==0  
-        saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_LRE',int2str(ifig)])
-        eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_LRE',int2str(ifig)]);
-        eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_LRE',int2str(ifig)]);
-        if options_.nograph, close(gcf); end
-    end
 end
 
-ifig=0;
-nbox = min(nparam-1,12);
-for j=1:nparam,  
-    if mod(j,12)==1,    
-        ifig = ifig+1;    
-        figure('name','Partial correlations in the model'),    
-        iplo=0;  
-    end
-    iplo=iplo+1;    
-    mmm = mean(squeeze(idemodel.Pco(:,j,:))');    
-    [sss, immm] = sort(-mmm);    
-    subplot(3,4,iplo), 
-    if nbox==1,
-        myboxplot(squeeze(idemodel.Pco(immm(2:nbox+1),j,:))),
-    else
-        myboxplot(squeeze(idemodel.Pco(immm(2:nbox+1),j,:))'),
-    end
-    set(gca,'ylim',[0 1])
-    set(gca,'xticklabel','')
-    for ip=1:nbox, %estim_params_.np,  
-        text(ip,-0.02,bayestopt_.name{immm(ip+1)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
-    end
-    title(bayestopt_.name{j}),
-    if j==nparam | mod(j,12)==0  
-        saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_model',int2str(ifig)])
-        eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_model',int2str(ifig)]);
-        eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_model',int2str(ifig)]);
-        if options_.nograph, close(gcf); end
-    end
-end
 
-ifig=0;
-nbox = min(nparam-1,12);
-for j=1:nparam,  
-    if mod(j,12)==1,    
-        ifig = ifig+1;    
-        figure('name','Partial correlations in the 1st and 2nd moments'),    
-        iplo=0;  
-    end
-    iplo=iplo+1;    
-    mmm = mean(squeeze(idemoments.Pco(:,j,:))');    
-    [sss, immm] = sort(-mmm);    
-    subplot(3,4,iplo), 
-    if nbox==1,
-        myboxplot(squeeze(idemoments.Pco(immm(2:nbox+1),j,:))),
+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
-        myboxplot(squeeze(idemoments.Pco(immm(2:nbox+1),j,:))'),
-    end
-    set(gca,'ylim',[0 1])
-    set(gca,'xticklabel','')
-    for ip=1:nbox, %estim_params_.np,  
-        text(ip,-0.02,bayestopt_.name{immm(ip+1)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
-    end
-    title(bayestopt_.name{j}),
-    if j==nparam | mod(j,12)==0  
-        saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_moments',int2str(ifig)])
-        eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_moments',int2str(ifig)]);
-        eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_moments',int2str(ifig)]);
-        if options_.nograph, close(gcf); end
+    pco(j+1:end,j) = abs(idelre.Pco(j+1:end,j))';    
     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
+% ifig=0;
+% nbox = min(np-1,12);
+% for j=1:np,  
+%     if mod(j,12)==1,    
+%         ifig = ifig+1;    
+%         figure('name','Pairwise correlations in the LRE model'),    
+%         iplo=0;  
+%     end
+%     iplo=iplo+1;    
+%     if SampleSize>1
+%     mmm = mean(abs(squeeze(idelre.Pco(:,j,:))'));    
+%     else
+%     mmm = abs(idelre.Pco(:,j))';    
+%     end
+%     [sss, immm] = sort(-mmm);    
+%     subplot(3,4,iplo), 
+%     if nbox==1,
+%         myboxplot(squeeze(idelre.Pco(immm(2:nbox+1),j,:))),
+%     else
+%         myboxplot(squeeze(idelre.Pco(immm(2:nbox+1),j,:))'),
+%     end
+%     set(gca,'ylim',[-1 1],'ygrid','on')
+%     set(gca,'xticklabel','')
+%     for ip=1:nbox, %np,  
+%         text(ip,-1.02,deblank(M_.param_names(indx(immm(ip+1)),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
+%     end
+%     title(deblank(M_.param_names(indx(j),:))),
+%     if j==np | mod(j,12)==0  
+%         saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_LRE',int2str(ifig)])
+%         eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_LRE',int2str(ifig)]);
+%         eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_LRE',int2str(ifig)]);
+%         if options_.nograph, close(gcf); end
+%     end
+% end
+% 
+% ifig=0;
+% nbox = min(nparam-1,12);
+% for j=1:nparam,  
+%     if mod(j,12)==1,    
+%         ifig = ifig+1;    
+%         figure('name','Pairwise correlations in the model'),    
+%         iplo=0;  
+%     end
+%     iplo=iplo+1;    
+%     if SampleSize>1,
+%     mmm = mean(abs(squeeze(idemodel.Pco(:,j,:))'));    
+%     else
+%     mmm = abs(idemodel.Pco(:,j)');    
+%     end
+%     [sss, immm] = sort(-mmm);    
+%     subplot(3,4,iplo), 
+%     if nbox==1,
+%         myboxplot(squeeze(idemodel.Pco(immm(2:nbox+1),j,:))),
+%     else
+%         myboxplot(squeeze(idemodel.Pco(immm(2:nbox+1),j,:))'),
+%     end
+%     set(gca,'ylim',[-1 1],'ygrid','on')
+%     set(gca,'xticklabel','')
+%     for ip=1:nbox, %np,  
+%         text(ip,-1.02,name{immm(ip+1)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+%     end
+%     title(name{j}),
+%     if j==nparam | mod(j,12)==0  
+%         saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_model',int2str(ifig)])
+%         eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_model',int2str(ifig)]);
+%         eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_model',int2str(ifig)]);
+%         if options_.nograph, close(gcf); end
+%     end
+% end
+% 
+% ifig=0;
+% nbox = min(nparam-1,12);
+% for j=1:nparam,  
+%     if mod(j,12)==1,    
+%         ifig = ifig+1;    
+%         figure('name','Pairwise correlations in the 1st and 2nd moments'),    
+%         iplo=0;  
+%     end
+%     iplo=iplo+1;    
+%     if SampleSize>1
+%     mmm = mean(abs(squeeze(idemoments.Pco(:,j,:))'));    
+%     else
+%     mmm = abs(idemoments.Pco(:,j)');    
+%     end
+%     [sss, immm] = sort(-mmm);    
+%     subplot(3,4,iplo), 
+%     if nbox==1,
+%         myboxplot(squeeze(idemoments.Pco(immm(2:nbox+1),j,:))),
+%     else
+%         myboxplot(squeeze(idemoments.Pco(immm(2:nbox+1),j,:))'),
+%     end
+%     set(gca,'ylim',[-1 1],'ygrid','on')
+%     set(gca,'xticklabel','')
+%     for ip=1:nbox, %np,  
+%         text(ip,-1.02,name{immm(ip+1)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+%     end
+%     title(name{j}),
+%     if j==nparam | mod(j,12)==0  
+%         saveas(gcf,[IdentifDirectoryName,'/',M_.fname,'_ident_PCORR_moments',int2str(ifig)])
+%         eval(['print -depsc2 ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_moments',int2str(ifig)]);
+%         eval(['print -dpdf ' IdentifDirectoryName '/' M_.fname '_ident_PCORR_moments',int2str(ifig)]);
+%         if options_.nograph, close(gcf); end
+%     end
+% end
 
 
 if exist('OCTAVE_VERSION')
@@ -625,3 +951,8 @@ if exist('OCTAVE_VERSION')
 else
     warning on MATLAB:dividebyzero
 end
+
+disp(' ')
+disp(['==== Identification analysis completed ====' ]),
+disp(' ')
+disp(' ')