diff --git a/matlab/identification_checks.m b/matlab/identification_checks.m
index 5da28b6b8..cce37622b 100644
--- a/matlab/identification_checks.m
+++ b/matlab/identification_checks.m
@@ -1,4 +1,4 @@
-function [McoH, McoJ, McoGP, PcoH, PcoJ, PcoGP, condH, condJ, condGP, eH, eJ, eGP, ind01, ind02, indnoH, indnoJ] = identification_checks(H,JJ, gp, bayestopt_)
+function [McoH, McoJ, McoGP, PcoH, PcoJ, PcoGP, condH, condJ, condGP, eH, eJ, eGP, ind01, ind02, indnoH, indnoJ, ixnoH, ixnoJ] = identification_checks(H,JJ, gp, bayestopt_)
 
 % Copyright (C) 2008 Dynare Team
 %
@@ -23,36 +23,42 @@ function [McoH, McoJ, McoGP, PcoH, PcoJ, PcoGP, condH, condJ, condGP, eH, eJ, eG
 % 1. check rank of H at theta
 npar = size(H,2);
 npar0 = size(gp,2);
-indnoH = {};
-indnoJ = {};
-indnoLRE = {};
-ind1 = find(vnorm(H)~=0);
+indnoH = zeros(1,npar);
+indnoJ = zeros(1,npar);
+indnoLRE = zeros(1,npar0);
+ind1 = find(vnorm(H)>=eps);
 H1 = H(:,ind1);
 covH = H1'*H1;
 sdH = sqrt(diag(covH));
 sdH = sdH*sdH';
-[e1,e2] = eig( (H1'*H1)./sdH );
+% [e1,e2] = eig( (H1'*H1)./sdH );
+[eu,e2,e1] = svd( H1, 0 );
 eH = zeros(npar,npar);
 % eH(ind1,:) = e1;
 eH(ind1,length(find(vnorm(H)==0))+1:end) = e1;
 eH(find(vnorm(H)==0),1:length(find(vnorm(H)==0)))=eye(length(find(vnorm(H)==0)));
 condH = cond(H1);
-% condH = cond(H1'*H1);
+condHH = cond(H1'*H1);
+rankH = rank(H);
+rankHH = rank(H'*H);
 
-ind2 = find(vnorm(JJ)~=0);
+ind2 = find(vnorm(JJ)>=eps);
 JJ1 = JJ(:,ind2);
 covJJ = JJ1'*JJ1;
 sdJJ = sqrt(diag(covJJ));
 sdJJ = sdJJ*sdJJ';
-[ee1,ee2] = eig( (JJ1'*JJ1)./sdJJ );
+% [ee1,ee2] = eig( (JJ1'*JJ1)./sdJJ );
+[eu,ee2,ee1] = svd( JJ1, 0 );
 % eJ = NaN(npar,length(ind2));
 eJ = zeros(npar,npar);
 eJ(ind2,length(find(vnorm(JJ)==0))+1:end) = ee1;
 eJ(find(vnorm(JJ)==0),1:length(find(vnorm(JJ)==0)))=eye(length(find(vnorm(JJ)==0)));
-% condJ = cond(JJ1'*JJ1);
+condJJ = cond(JJ1'*JJ1);
 condJ = cond(JJ1);
+rankJJ = rank(JJ'*JJ);
+rankJ = rank(JJ);
 
-ind3 = find(vnorm(gp)~=0);
+ind3 = find(vnorm(gp)>=eps);
 gp1 = gp(:,ind3);
 covgp = gp1'*gp1;
 sdgp = sqrt(diag(covgp));
@@ -65,54 +71,11 @@ eGP(find(vnorm(gp)==0),1:length(find(vnorm(gp)==0)))=eye(length(find(vnorm(gp)==
 % condJ = cond(JJ1'*JJ1);
 condGP = cond(gp1);
 
-if rank(H)<npar
-    ixno = 0;
-    %         - find out which parameters are involved,
-    % using something like the vnorm and the eigenvalue decomposition of H;
-    %   disp('Some parameters are NOT identified in the model: H rank deficient')
-    %   disp(' ')
-    if length(ind1)<npar,
-        ixno = ixno + 1;
-        indnoH(ixno) = {find(~ismember([1:npar],ind1))};
-        %     disp('Not identified params')
-        %     disp(bayestopt_.name(indnoH{1}))
-        %     disp(' ')
-    end
-    e0 = find(abs(diag(e2))<eps);
-    for j=1:length(e0),
-        ixno = ixno + 1;
-        indnoH(ixno) = {ind1(find(e1(:,e0(j))))};
-        %     disp('Perfectly collinear parameters')
-        %     disp(bayestopt_.name(indnoH{ixno}))
-        %     disp(' ')
-    end
-else % rank(H)==length(theta), go to 2
-     % 2. check rank of J
-     %   disp('All parameters are identified at theta in the model (rank of H)')
-     %   disp(' ')
-end
-
-if rank(JJ)<npar
-    ixno = 0;
-    %         - find out which parameters are involved
-    %   disp('Some parameters are NOT identified by the moments included in J')
-    %   disp(' ')
-    if length(ind2)<npar,
-        ixno = ixno + 1;
-        indnoJ(ixno) = {find(~ismember([1:npar],ind2))};
-    end
-    ee0 = find(abs(diag(ee2))<eps);
-    for j=1:length(ee0),
-        ixno = ixno + 1;
-        indnoJ(ixno) = {ind2(find(ee1(:,ee0(j))))};
-        %     disp('Perfectly collinear parameters in moments J')
-        %     disp(bayestopt_.name(indnoJ{ixno}))
-        %     disp(' ')
-    end
-else  %rank(J)==length(theta) =>
-      %   disp('All parameters are identified at theta by the moments included in J')
-end
 
+ind01 = zeros(npar,1);
+ind02 = zeros(npar,1);
+ind01(ind1) = 1;
+ind02(ind2) = 1;
 
 % rank(H1)==size(H1,2)
 % rank(JJ1)==size(JJ1,2)
@@ -135,6 +98,64 @@ end
 % format long  % some are nearly 1
 % McoJ
 
+ixno = 0;
+if rankH<npar | rankHH<npar | min(1-McoH)<1.e-10
+    %         - find out which parameters are involved,
+    % using something like the vnorm and the eigenvalue decomposition of H;
+    %   disp('Some parameters are NOT identified in the model: H rank deficient')
+    %   disp(' ')
+    if length(ind1)<npar,
+        ixno = ixno + 1;
+%         indnoH(ixno) = {find(~ismember([1:npar],ind1))};
+        indnoH(ixno,:) = (~ismember([1:npar],ind1));
+        %     disp('Not identified params')
+        %     disp(bayestopt_.name(indnoH{1}))
+        %     disp(' ')
+    end
+    e0 = [rankHH+1:length(ind1)];
+    for j=1:length(e0),
+        ixno = ixno + 1;
+%         indnoH(ixno) = {ind1(find(abs(e1(:,e0(j)))) > 1.e-6 )};
+        indnoH(ixno,:) = (abs(e1(:,e0(j))) > 1.e-6 )';
+        %     disp('Perfectly collinear parameters')
+        %     disp(bayestopt_.name(indnoH{ixno}))
+        %     disp(' ')
+%         ind01(indnoH{ixno})=0;
+    end
+else % rank(H)==length(theta), go to 2
+     % 2. check rank of J
+     %   disp('All parameters are identified at theta in the model (rank of H)')
+     %   disp(' ')
+end
+ixnoH=ixno;
+
+ixno = 0;
+if rankJ<npar | rankJJ<npar | min(1-McoJ)<1.e-10
+    %         - find out which parameters are involved
+    %   disp('Some parameters are NOT identified by the moments included in J')
+    %   disp(' ')
+    if length(ind2)<npar,
+        ixno = ixno + 1;
+%         indnoJ(ixno) = {find(~ismember([1:npar],ind2))};
+        indnoJ(ixno,:) = (~ismember([1:npar],ind2));
+    end
+    ee0 = [rankJJ+1:length(ind2)];
+    if isempty(ee0),
+        cccc=0';
+    end
+    for j=1:length(ee0),
+        ixno = ixno + 1;
+%         indnoJ(ixno) = {ind2( find(abs(ee1(:,ee0(j))) > 1.e-6) )};
+        indnoJ(ixno,:) = (abs(ee1(:,ee0(j))) > 1.e-6)';
+        %     disp('Perfectly collinear parameters in moments J')
+        %     disp(bayestopt_.name(indnoJ{ixno}))
+        %     disp(' ')
+%         ind02(indnoJ{ixno})=0;
+    end
+else  %rank(J)==length(theta) =>
+%         disp('All parameters are identified at theta by the moments included in J')
+end
+ixnoJ=ixno;
 
 % here there is no exact linear dependence, but there are several
 %     near-dependencies, mostly due to strong pairwise colliniearities, which can
@@ -168,10 +189,10 @@ for ii = 1:size(gp1,2);
 end
 
 
-ind01 = zeros(npar,1);
-ind02 = zeros(npar,1);
-ind01(ind1) = 1;
-ind02(ind2) = 1;
+% ind01 = zeros(npar,1);
+% ind02 = zeros(npar,1);
+% ind01(ind1) = 1;
+% ind02(ind2) = 1;