diff --git a/matlab/disp_identification.m b/matlab/disp_identification.m
index 25cbb114e797a878d2546faf43401fd4e0d974cd..1c40f28094980b3b5cf85dfdf5a7235c5d2d1010 100644
--- a/matlab/disp_identification.m
+++ b/matlab/disp_identification.m
@@ -24,7 +24,7 @@ function disp_identification(pdraws, ide_reducedform, ide_moments, ide_spectrum,
 % This function is called by
 %   * dynare_identification.m
 % =========================================================================
-% Copyright (C) 2010-2019 Dynare Team
+% Copyright (C) 2010-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -82,7 +82,7 @@ for jide = 1:4
     %% Set output strings depending on test
     if jide == 1
         strTest = 'REDUCED-FORM'; strJacobian = 'Tau'; strMeaning = 'Jacobian of steady state and reduced-form solution matrices';
-        if ~no_identification_reducedform
+        if ~no_identification_reducedform && ~ isempty(fieldnames(ide_reducedform))
             noidentification = 0; ide = ide_reducedform;
             if SampleSize == 1
                 Jacob = ide.dREDUCEDFORM;
@@ -97,7 +97,7 @@ for jide = 1:4
         elseif options_ident.order == 3
             strMeaning = 'Jacobian of first-order minimal system and third-order accurate mean';
         end
-        if ~no_identification_minimal
+        if ~no_identification_minimal && ~(length(ide_minimal.minimal_state_space)==1 && ide_minimal.minimal_state_space==0)
             noidentification = 0; ide = ide_minimal;
             if SampleSize == 1
                 Jacob = ide.dMINIMAL;
@@ -110,7 +110,7 @@ for jide = 1:4
         if options_ident.order > 1
             strTest = 'SPECTRUM (Mutschler, 2015)';
         end
-        if ~no_identification_spectrum
+        if ~no_identification_spectrum && ~isempty(fieldnames(ide_spectrum))
             noidentification = 0; ide = ide_spectrum;
             if SampleSize == 1
                 Jacob = ide.dSPECTRUM;
@@ -123,7 +123,7 @@ for jide = 1:4
         if options_ident.order > 1
             strTest = 'MOMENTS (Mutschler, 2015)'; strJacobian = 'Mbar';
         end
-        if ~no_identification_moments
+        if ~no_identification_moments && ~isempty(fieldnames(ide_moments))
             noidentification = 0; ide = ide_moments;
             if SampleSize == 1
                 Jacob = ide.si_dMOMENTS;
@@ -136,20 +136,39 @@ for jide = 1:4
     if ~noidentification
         %% display problematic parameters computed by identifcation_checks.m
         if ~checks_via_subsets
-            if any(ide.ino) || any(any(ide.ind0==0)) || any(any(ide.jweak_pair))
+            EffectiveSampleSize=SampleSize;
+            non_minimal_state_space_error=0;
+            if SampleSize>1 && jide==2  && any(~ide.minimal_state_space)
+                EffectiveSampleSize=SampleSize-sum(~ide.minimal_state_space);
+                non_minimal_state_space_error=1;
+                ide.ino(~ide.minimal_state_space,:)=[];
+                ide.ind0(~ide.minimal_state_space,:)=[];
+                ide.jweak_pair(~ide.minimal_state_space,:)=[];
+                ide.minimal_state_space(~ide.minimal_state_space,:)=[];
+            end
+            if any(ide.ino) || any(any(ide.ind0==0)) || any(any(ide.jweak_pair)) || non_minimal_state_space_error
                 no_warning_message_display=0;
                 skipline()
                 disp([upper(strTest), ':'])
                 disp('    !!!WARNING!!!');
                 if SampleSize>1
-                    disp(['    The rank of ', strJacobian, ' (', strMeaning, ') is deficient for ', num2str(length(find(ide.ino))),' out of ',int2str(SampleSize),' MC runs!'  ]),
+                    if non_minimal_state_space_error
+                        fprintf(['\n    The minimal state space could not be computed for %u out of %u cases.\n'],SampleSize-EffectiveSampleSize,SampleSize);
+                    end
+                    if jide==2
+                        if sum(ide.ino & ide.minimal_state_space)>0
+                        disp(['    The rank of ', strJacobian, ' (', strMeaning, ') is deficient for ', num2str(sum(ide.ino & ide.minimal_state_space)),' out of ',int2str(EffectiveSampleSize),' effective MC runs!'  ])
+                        end
+                    else
+                        disp(['    The rank of ', strJacobian, ' (', strMeaning, ') is deficient for ', num2str(sum(ide.ino)),' out of ',int2str(EffectiveSampleSize),' effective MC runs!'  ]),
+                    end
                 else
                     disp(['    The rank of ', strJacobian, ' (', strMeaning, ') is deficient!']),
                 end
                 skipline()
                 for j=1:totparam_nbr
                     if any(ide.ind0(:,j)==0)
-                        pno = 100*length(find(ide.ind0(:,j)==0))/SampleSize;
+                        pno = 100*length(find(ide.ind0(:,j)==0))/EffectiveSampleSize;
                         if SampleSize>1
                             disp(['    ',name{j},' is not identified for ',num2str(pno),'% of MC runs!' ])
                         else
@@ -166,7 +185,7 @@ for jide = 1:4
                         [jx,jy]=find(jmap_pair==j);
                         jstore=[jstore jx(1) jy(1)];
                         if SampleSize > 1
-                            disp(['    [',name{jx(1)},',',name{jy(1)},'] are PAIRWISE collinear for ',num2str((iweak)/SampleSize*100),'% of MC runs!' ])
+                            disp(['    [',name{jx(1)},',',name{jy(1)},'] are PAIRWISE collinear for ',num2str((iweak)/EffectiveSampleSize*100),'% of MC runs!' ])
                         else
                             disp(['    [',name{jx(1)},',',name{jy(1)},'] are PAIRWISE collinear!' ])
                         end
@@ -176,7 +195,7 @@ for jide = 1:4
                     iweak = length(find(ide.jweak(:,j)));
                     if iweak && ~ismember(j,jstore)
                         if SampleSize>1
-                            disp(['    ',name{j},' is collinear w.r.t. all other parameters for ',num2str(iweak/SampleSize*100),'% of MC runs!' ])
+                            disp(['    ',name{j},' is collinear w.r.t. all other parameters for ',num2str(iweak/EffectiveSampleSize*100),'% of MC runs!' ])
                         else
                             disp(['    ',name{j},' is collinear w.r.t. all other parameters!' ])
                         end
@@ -241,7 +260,7 @@ end
 
 
 
-%% Advanced identificaton patterns
+%% Advanced identification patterns
 if SampleSize==1 && options_ident.advanced
     skipline()
     for j=1:size(ide_moments.cosndMOMENTS,2)
diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m
index a631600e5f651027469749e3154a78e2b6a90a55..ab39a7242e505235637ed80c8b9832a98ff6e913 100644
--- a/matlab/dynare_identification.m
+++ b/matlab/dynare_identification.m
@@ -481,7 +481,7 @@ if iload <=0
     end
     options_ident.tittxt = parameters; %title text for graphs and figures
     % perform identification analysis for single point
-    [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, options_ident] = ...
+    [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info] = ...
         identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end implies initialization of persistent variables
     if info(1)~=0
         % there are errors in the solution algorithm
@@ -498,7 +498,7 @@ if iload <=0
                 params = prior_draw();
                 options_ident.tittxt = 'Random_prior_params'; %title text for graphs and figures
                 % perform identification analysis
-                [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, options_ident] = ...
+                [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info] = ...
                     identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1);
             end
         end
@@ -552,7 +552,7 @@ if iload <=0
         end
         options_ident.tittxt = []; % clear title text for graphs and figures
         % run identification analysis
-        [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, ide_derivatives_info, info, options_MC] = ...
+        [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, ide_derivatives_info, info, error_indicator] = ...
             identification_analysis(params, indpmodel, indpstderr, indpcorr, options_MC, dataset_info, prior_exist, 0); % the 0 implies that we do not initialize persistent variables anymore
 
         if iteration==0 && info(1)==0 % preallocate storage in the first admissable run
@@ -564,7 +564,7 @@ if iload <=0
             STO_si_dDYNAMIC          = zeros([size(ide_dynamic.si_dDYNAMIC, 1), modparam_nbr, MAX_RUNS_BEFORE_SAVE_TO_FILE]);
             STO_DYNAMIC              = zeros(size(ide_dynamic.DYNAMIC, 1), SampleSize);
             IDE_DYNAMIC.ind_dDYNAMIC = ide_dynamic.ind_dDYNAMIC;
-            IDE_DYNAMIC.in0          = zeros(SampleSize, modparam_nbr);
+            IDE_DYNAMIC.ino          = zeros(SampleSize, modparam_nbr);
             IDE_DYNAMIC.ind0         = zeros(SampleSize, modparam_nbr);
             IDE_DYNAMIC.jweak        = zeros(SampleSize, modparam_nbr);
             IDE_DYNAMIC.jweak_pair   = zeros(SampleSize, modparam_nbr*(modparam_nbr+1)/2);
@@ -576,7 +576,7 @@ if iload <=0
                 STO_si_dREDUCEDFORM              = zeros([size(ide_reducedform.si_dREDUCEDFORM, 1), totparam_nbr, MAX_RUNS_BEFORE_SAVE_TO_FILE]);
                 STO_REDUCEDFORM                  = zeros(size(ide_reducedform.REDUCEDFORM, 1), SampleSize);
                 IDE_REDUCEDFORM.ind_dREDUCEDFORM = ide_reducedform.ind_dREDUCEDFORM;
-                IDE_REDUCEDFORM.in0              = zeros(SampleSize, 1);
+                IDE_REDUCEDFORM.ino              = zeros(SampleSize, 1);
                 IDE_REDUCEDFORM.ind0             = zeros(SampleSize, totparam_nbr);
                 IDE_REDUCEDFORM.jweak            = zeros(SampleSize, totparam_nbr);
                 IDE_REDUCEDFORM.jweak_pair       = zeros(SampleSize, totparam_nbr*(totparam_nbr+1)/2);
@@ -593,7 +593,7 @@ if iload <=0
                 STO_si_dMOMENTS          = zeros([size(ide_moments.si_dMOMENTS, 1), totparam_nbr, MAX_RUNS_BEFORE_SAVE_TO_FILE]);
                 STO_MOMENTS              = zeros(size(ide_moments.MOMENTS, 1), SampleSize);
                 IDE_MOMENTS.ind_dMOMENTS = ide_moments.ind_dMOMENTS;
-                IDE_MOMENTS.in0          = zeros(SampleSize, 1);
+                IDE_MOMENTS.ino          = zeros(SampleSize, 1);
                 IDE_MOMENTS.ind0         = zeros(SampleSize, totparam_nbr);
                 IDE_MOMENTS.jweak        = zeros(SampleSize, totparam_nbr);
                 IDE_MOMENTS.jweak_pair   = zeros(SampleSize, totparam_nbr*(totparam_nbr+1)/2);
@@ -611,7 +611,7 @@ if iload <=0
             if ~options_MC.no_identification_spectrum
                 STO_dSPECTRUM              = zeros([size(ide_spectrum.dSPECTRUM, 1), size(ide_spectrum.dSPECTRUM, 2), MAX_RUNS_BEFORE_SAVE_TO_FILE]);
                 IDE_SPECTRUM.ind_dSPECTRUM = ide_spectrum.ind_dSPECTRUM;
-                IDE_SPECTRUM.in0           = zeros(SampleSize, 1);
+                IDE_SPECTRUM.ino           = zeros(SampleSize, 1);
                 IDE_SPECTRUM.ind0          = zeros(SampleSize, totparam_nbr);
                 IDE_SPECTRUM.jweak         = zeros(SampleSize, totparam_nbr);
                 IDE_SPECTRUM.jweak_pair    = zeros(SampleSize, totparam_nbr*(totparam_nbr+1)/2);
@@ -626,12 +626,13 @@ if iload <=0
             if ~options_MC.no_identification_minimal
                 STO_dMINIMAL             = zeros([size(ide_minimal.dMINIMAL, 1), size(ide_minimal.dMINIMAL, 2), MAX_RUNS_BEFORE_SAVE_TO_FILE]);
                 IDE_MINIMAL.ind_dMINIMAL = ide_minimal.ind_dMINIMAL;
-                IDE_MINIMAL.in0          = zeros(SampleSize, 1);
+                IDE_MINIMAL.ino          = zeros(SampleSize, 1);
                 IDE_MINIMAL.ind0         = zeros(SampleSize, totparam_nbr);
                 IDE_MINIMAL.jweak        = zeros(SampleSize, totparam_nbr);
                 IDE_MINIMAL.jweak_pair   = zeros(SampleSize, totparam_nbr*(totparam_nbr+1)/2);
                 IDE_MINIMAL.cond         = zeros(SampleSize, 1);
                 IDE_MINIMAL.Mco          = zeros(SampleSize, totparam_nbr);
+                IDE_MINIMAL.minimal_state_space = zeros(SampleSize, 1);
             else
                 STO_dMINIMAL = {};
                 IDE_MINIMAL  = {};
@@ -654,7 +655,7 @@ if iload <=0
             IDE_DYNAMIC.Mco(iteration,:)        = ide_dynamic.Mco;
 
             % store results for reduced form
-            if ~options_MC.no_identification_reducedform
+            if ~options_MC.no_identification_reducedform && ~error_indicator.identification_reducedform
                 STO_REDUCEDFORM(:,iteration)            = ide_reducedform.REDUCEDFORM;
                 STO_si_dREDUCEDFORM(:,:,run_index)      = ide_reducedform.si_dREDUCEDFORM;
                 IDE_REDUCEDFORM.cond(iteration,1)       = ide_reducedform.cond;
@@ -666,7 +667,7 @@ if iload <=0
             end
 
             % store results for moments
-            if ~options_MC.no_identification_moments
+            if ~options_MC.no_identification_moments && ~error_indicator.identification_moments
                 STO_MOMENTS(:,iteration)            = ide_moments.MOMENTS;
                 STO_si_dMOMENTS(:,:,run_index)      = ide_moments.si_dMOMENTS;
                 IDE_MOMENTS.cond(iteration,1)       = ide_moments.cond;
@@ -680,7 +681,7 @@ if iload <=0
             end
 
             % store results for spectrum
-            if ~options_MC.no_identification_spectrum
+            if ~options_MC.no_identification_spectrum && ~error_indicator.identification_spectrum
                 STO_dSPECTRUM(:,:,run_index)         = ide_spectrum.dSPECTRUM;
                 IDE_SPECTRUM.cond(iteration,1)       = ide_spectrum.cond;
                 IDE_SPECTRUM.ino(iteration,1)        = ide_spectrum.ino;
@@ -691,14 +692,19 @@ if iload <=0
             end
 
             % store results for minimal system
-            if ~options_MC.no_identification_minimal
-                STO_dMINIMAL(:,:,run_index)         = ide_minimal.dMINIMAL;
-                IDE_MINIMAL.cond(iteration,1)       = ide_minimal.cond;
-                IDE_MINIMAL.ino(iteration,1)        = ide_minimal.ino;
-                IDE_MINIMAL.ind0(iteration,:)       = ide_minimal.ind0;
-                IDE_MINIMAL.jweak(iteration,:)      = ide_minimal.jweak;
-                IDE_MINIMAL.jweak_pair(iteration,:) = ide_minimal.jweak_pair;
-                IDE_MINIMAL.Mco(iteration,:)        = ide_minimal.Mco;
+            if ~options_MC.no_identification_minimal 
+                if ~error_indicator.identification_minimal
+                    STO_dMINIMAL(:,:,run_index)                  = ide_minimal.dMINIMAL;
+                    IDE_MINIMAL.cond(iteration,1)                = ide_minimal.cond;
+                    IDE_MINIMAL.ino(iteration,1)                 = ide_minimal.ino;
+                    IDE_MINIMAL.ind0(iteration,:)                = ide_minimal.ind0;
+                    IDE_MINIMAL.jweak(iteration,:)               = ide_minimal.jweak;
+                    IDE_MINIMAL.jweak_pair(iteration,:)          = ide_minimal.jweak_pair;
+                    IDE_MINIMAL.Mco(iteration,:)                 = ide_minimal.Mco;
+                    IDE_MINIMAL.minimal_state_space(iteration,:) = ide_minimal.minimal_state_space;
+                else
+                    IDE_MINIMAL.minimal_state_space(iteration,:) = ide_minimal.minimal_state_space;
+                end
             end
 
             % save results to file: either to avoid running into memory issues, i.e. (run_index==MAX_RUNS_BEFORE_SAVE_TO_FILE) or if finished (iteration==SampleSize)
@@ -900,7 +906,7 @@ if SampleSize > 1
                 fprintf('\nTesting %s.\n',tittxt);
                 if ~iload
                     options_ident.tittxt = tittxt; %title text for graphs and figures
-                    [ide_moments_max, ide_spectrum_max, ide_minimal_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, derivatives_info_max, info_max, options_ident] = ...
+                    [ide_moments_max, ide_spectrum_max, ide_minimal_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, derivatives_info_max, info_max] = ...
                         identification_analysis(pdraws(jmax,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes some persistent variables
                     save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_max', 'ide_moments_max', 'ide_spectrum_max', 'ide_minimal_max','ide_reducedform_max', 'ide_dynamic_max', 'jmax', '-append');
                 end
@@ -918,7 +924,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, options_ident] = ...
+                    [ide_moments_min, ide_spectrum_min, ide_minimal_min, ide_hess_min, ide_reducedform_min, ide_dynamic_min, derivatives_info_min, info_min] = ...
                         identification_analysis(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
@@ -940,7 +946,7 @@ if SampleSize > 1
                     fprintf('\nTesting %s.\n',tittxt);
                     if ~iload
                         options_ident.tittxt = tittxt; %title text for graphs and figures
-                        [ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), derivatives_info_(j), info_resolve, options_ident] = ...
+                        [ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), derivatives_info_(j), info_resolve] = ...
                             identification_analysis(pdraws(jcrit(j),:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1);
                     end
                     advanced0 = options_ident.advanced; options_ident.advanced = 1; %make sure advanced setting is on
diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m
index fd6c2fef4039679b25b8160d8e5ab85a0d189fa4..70af85aa22ffee84ca007547d796491d31128b69 100644
--- a/matlab/identification_analysis.m
+++ b/matlab/identification_analysis.m
@@ -1,5 +1,5 @@
-function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, options_ident] = identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init)
-% [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, options_ident] = identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init)
+function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, error_indicator] = identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init)
+% [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, error_indicator] = identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init)
 % -------------------------------------------------------------------------
 % This function wraps all identification analysis, i.e. it
 % (1) wraps functions for the theoretical identification analysis based on moments (Iskrev, 2010),
@@ -47,8 +47,8 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
 %                         info about first-order perturbation derivatives, used in dsge_likelihood.m
 %    * info               [integer]
 %                         output from dynare_resolve
-%    * options_ident      [structure]
-%                         updated identification options
+%    * error_indicator    [structure]
+%                         indicator on problems
 % -------------------------------------------------------------------------
 % This function is called by
 %   * dynare_identification.m
@@ -71,7 +71,7 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
 %   * stoch_simul
 %   * vec
 % =========================================================================
-% Copyright (C) 2008-2020 Dynare Team
+% Copyright (C) 2008-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -115,7 +115,6 @@ if ~isempty(estim_params_)
 end
 
 %get options (see dynare_identification.m for description of options)
-order               = options_ident.order;
 nlags               = options_ident.ar;
 advanced            = options_ident.advanced;
 replic              = options_ident.replic;
@@ -126,11 +125,11 @@ checks_via_subsets  = options_ident.checks_via_subsets;
 tol_deriv           = options_ident.tol_deriv;
 tol_rank            = options_ident.tol_rank;
 tol_sv              = options_ident.tol_sv;
-no_identification_strength    = options_ident.no_identification_strength;
-no_identification_reducedform = options_ident.no_identification_reducedform;
-no_identification_moments     = options_ident.no_identification_moments;
-no_identification_minimal     = options_ident.no_identification_minimal;
-no_identification_spectrum    = options_ident.no_identification_spectrum;
+error_indicator.identification_strength=0;
+error_indicator.identification_reducedform=0;
+error_indicator.identification_moments=0;
+error_indicator.identification_minimal=0;
+error_indicator.identification_spectrum=0;
 
 %Compute linear approximation and fill dr structure
 [oo_.dr,info,M_,oo_] = compute_decision_rules(M_,options_,oo_);
@@ -140,71 +139,71 @@ if info(1) == 0 %no errors in solution
     [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params_, M_, oo_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs);
     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
-        no_identification_minimal = 1;
-        options_ident.no_identification_minimal = 1;
+        error_indicator.identification_minimal = 1;
+        %options_ident.no_identification_minimal = 1;
     end
 
     if init
         %check stationarity
-        if ~no_identification_moments
+        if ~options_ident.no_identification_moments
             if any(any(isnan(MOMENTS)))
                 if options_.diffuse_filter == 1 % use options_ as it inherits diffuse_filter from options_ident if set by user
                     error('There are NaN''s in the theoretical moments. Make sure that for non-stationary models stationary transformations of non-stationary observables are used for checking identification. [TIP: use first differences].')
                 else
                     error('There are NaN''s in the theoretical moments. Please check whether your model has units roots, and you forgot to set diffuse_filter=1.' )
                 end
+                error_indicator.identification_moments=1;
             end
             ind_dMOMENTS = (find(max(abs(dMOMENTS'),[],1) > tol_deriv)); %index for non-zero rows
             if isempty(ind_dMOMENTS) && any(any(isnan(dMOMENTS)))                
-                error('There are NaN in the dMOMENTS matrix.' )
-            end
-            
+                error('There are NaN in the dMOMENTS matrix.')
+                error_indicator.identification_moments=1;
+            end            
         end
-        if ~no_identification_spectrum
+        if ~options_ident.no_identification_spectrum
             ind_dSPECTRUM = (find(max(abs(dSPECTRUM'),[],1) > tol_deriv)); %index for non-zero rows
             if any(any(isnan(dSPECTRUM)))
                 warning_SPECTRUM = 'WARNING: There are NaN in the dSPECTRUM matrix. Note that identification based on spectrum does not support non-stationary models (yet).\n';
                 warning_SPECTRUM = [warning_SPECTRUM '         Skip identification analysis based on spectrum.\n'];
                 fprintf(warning_SPECTRUM);
-                %reset options to neither display nor plot dSPECTRUM anymore
-                no_identification_spectrum = 1;
-                options_ident.no_identification_spectrum = 1;
+                %set indicator to neither display nor plot dSPECTRUM anymore
+                error_indicator.identification_spectrum = 1;                
             end
         end
-        if ~no_identification_minimal
+        if ~options_ident.no_identification_minimal
             ind_dMINIMAL = (find(max(abs(dMINIMAL'),[],1) > tol_deriv)); %index for non-zero rows
             if any(any(isnan(dMINIMAL)))
                 warning_MINIMAL = 'WARNING: There are NaN in the dMINIMAL matrix. Note that identification based on minimal system does not support non-stationary models (yet).\n';
                 warning_MINIMAL = [warning_MINIMAL '         Skip identification analysis based on minimal system.\n'];
                 fprintf(warning_MINIMAL);
-                %reset options to neither display nor plot dMINIMAL anymore
-                no_identification_minimal = 1;
-                options_ident.no_identification_minimal = 1;
+                %set indicator to neither display nor plot dMINIMAL anymore
+                error_indicator.identification_minimal = 1;
             end
         end
-        if no_identification_moments && no_identification_minimal && no_identification_spectrum
+        %The following cannot be reached yet due to erroring out when
+        %error_indicator.identification_moments is triggered
+        if error_indicator.identification_moments && error_indicator.identification_minimal && error_indicator.identification_spectrum
             %display error if all three criteria fail            
             error(sprintf('identification_analyis: Stationarity condition(s) failed and/or diffuse_filter option missing.\nMake sure that for non-stationary models stationary transformations of non-stationary observables are used for checking identification.\n[TIP: use first differences].'));
         end
 
         % Check order conditions
-        if ~no_identification_moments
+        if ~options_ident.no_identification_moments && ~error_indicator.identification_moments
             %check order condition of Iskrev (2010)
             while length(ind_dMOMENTS) < totparam_nbr && nlags < 10
                 %Try to add lags to autocovariogram if order condition fails
                 disp('The number of moments with non-zero derivative is smaller than the number of parameters')
                 disp(['Try increasing ar = ', int2str(nlags+1)])
                 nlags = nlags + 1;
-                options_ident.no_identification_minimal  = 1; %do not recompute dMINIMAL
-                options_ident.no_identification_spectrum = 1; %do not recompute dSPECTRUM
-                options_ident.ar = nlags;     %store new lag number
-                options_.ar      = nlags;     %store new lag number
-                [~, ~, ~, ~, ~, ~, MOMENTS, dMOMENTS, ~, ~, ~, ~] = get_identification_jacobians(estim_params_, M_, oo_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs);
+                options_ident_local=options_ident;
+                options_ident_local.no_identification_minimal  = 1; %do not recompute dMINIMAL
+                options_ident_local.no_identification_spectrum = 1; %do not recompute dSPECTRUM
+                options_ident_local.ar = nlags;     %store new lag number
+                options_.ar      = nlags;           %store new lag number
+                [~, ~, ~, ~, ~, ~, MOMENTS, dMOMENTS, ~, ~, ~, ~] = get_identification_jacobians(estim_params_, M_, oo_, options_, options_ident_local, indpmodel, indpstderr, indpcorr, indvobs);
 
                 ind_dMOMENTS = (find(max(abs(dMOMENTS'),[],1) > tol_deriv)); %new index with non-zero rows
             end
-            options_ident.no_identification_minimal  = no_identification_minimal;  % reset option to original setting
-            options_ident.no_identification_spectrum = no_identification_spectrum; % reset option to original setting
             if length(ind_dMOMENTS) < totparam_nbr
                 warning_MOMENTS = 'WARNING: Order condition for dMOMENTS failed: There are not enough moments and too many parameters.\n';
                 warning_MOMENTS = [warning_MOMENTS '         The number of moments with non-zero derivative is smaller than the number of parameters up to 10 lags.\n'];
@@ -212,31 +211,30 @@ if info(1) == 0 %no errors in solution
                 warning_MOMENTS = [warning_MOMENTS '         Skip identification analysis based on moments.\n'];
                 warning_MOMENTS = [warning_MOMENTS '         Skip identification strenght analysis.\n'];
                 fprintf(warning_MOMENTS);
-                %reset options to neither display nor plot dMOMENTS anymore
-                no_identification_moments = 1;
-                options_ident.no_identification_moments = 1;
-                no_identification_strength = 1;
-                options_ident.no_identification_strength = 1;
+                %set indicator to neither display nor plot dMOMENTS anymore
+                error_indicator.identification_moments = 1;
+                %options_ident.no_identification_moments = 1;
+                error_indicator.identification_strength = 1;
+                %options_ident.no_identification_strength = 1;
             end
         end
-        if ~no_identification_minimal
+        if ~options_ident.no_identification_minimal && ~error_indicator.identification_minimal
             if length(ind_dMINIMAL) < size(dMINIMAL,2)
                 warning_MINIMAL = 'WARNING: Order condition for dMINIMAL failed: There are too many parameters or too few observable variables.\n';
                 warning_MINIMAL = [warning_MINIMAL '         The number of minimal system elements with non-zero derivative is smaller than the number of parameters.\n'];
                 warning_MINIMAL = [warning_MINIMAL '         Either reduce the list of parameters, or increase number of observables.\n'];
                 warning_MINIMAL = [warning_MINIMAL '         Skip identification analysis based on minimal state space system.\n'];
                 fprintf(warning_MINIMAL);
-                %resest options to neither display nor plot dMINIMAL anymore
-                no_identification_minimal = 1;
-                options_ident.no_identification_minimal = 1;
+                %set indicator to neither display nor plot dMINIMAL anymore
+                error_indicator.identification_minimal = 1;                
             end
         end
         %Note that there is no order condition for dSPECTRUM, as the matrix is always of dimension totparam_nbr by totparam_nbr
-        if no_identification_moments && no_identification_minimal && no_identification_spectrum
+        if error_indicator.identification_moments && error_indicator.identification_minimal && error_indicator.identification_spectrum
             %error if all three criteria fail
             error('identification_analyis: Order condition(s) failed');
         end
-        if ~no_identification_reducedform
+        if ~options_ident.no_identification_reducedform && ~error_indicator.identification_reducedform
             ind_dREDUCEDFORM = (find(max(abs(dREDUCEDFORM'),[],1) > tol_deriv)); %index with non-zero rows
         end
         ind_dDYNAMIC = (find(max(abs(dDYNAMIC'),[],1) > tol_deriv)); %index with non-zero rows
@@ -244,16 +242,16 @@ if info(1) == 0 %no errors in solution
 
     DYNAMIC = DYNAMIC(ind_dDYNAMIC); %focus only on non-zero entries
     si_dDYNAMIC = (dDYNAMIC(ind_dDYNAMIC,:)); %focus only on non-zero rows
-    if ~no_identification_reducedform
+    if ~options_ident.no_identification_reducedform && ~error_indicator.identification_reducedform
         REDUCEDFORM = REDUCEDFORM(ind_dREDUCEDFORM); %focus only on non-zero entries
         si_dREDUCEDFORM = (dREDUCEDFORM(ind_dREDUCEDFORM,:)); %focus only on non-zero rows
     end
 
-    if ~no_identification_moments
+    if ~options_ident.no_identification_moments && ~error_indicator.identification_moments
         MOMENTS = MOMENTS(ind_dMOMENTS); %focus only on non-zero entries
         si_dMOMENTS   = (dMOMENTS(ind_dMOMENTS,:)); %focus only on non-zero derivatives
-%% MOMENTS IDENTIFICATION STRENGTH ANALYSIS
-        if ~no_identification_strength && init %only for initialization of persistent vars
+        %% MOMENTS IDENTIFICATION STRENGTH ANALYSIS
+        if ~options_ident.no_identification_strength && ~error_indicator.identification_strength && init %only for initialization of persistent vars
             ide_strength_dMOMENTS        = NaN(1,totparam_nbr); %initialize
             ide_strength_dMOMENTS_prior  = NaN(1,totparam_nbr); %initialize
             ide_uncert_unnormaliz = NaN(1,totparam_nbr); %initialize
@@ -438,7 +436,7 @@ if info(1) == 0 %no errors in solution
     ide_dynamic.dDYNAMIC      = dDYNAMIC;
     ide_dynamic.DYNAMIC       = DYNAMIC;
 
-    if ~no_identification_reducedform
+    if ~options_ident.no_identification_reducedform && ~error_indicator.identification_reducedform
         if normalize_jacobians
             norm_dREDUCEDFORM = max(abs(si_dREDUCEDFORM),[],2);
             norm_dREDUCEDFORM = norm_dREDUCEDFORM(:,ones(totparam_nbr,1));
@@ -453,7 +451,7 @@ if info(1) == 0 %no errors in solution
         ide_reducedform.REDUCEDFORM       = REDUCEDFORM;
     end
 
-    if ~no_identification_moments
+    if ~options_ident.no_identification_moments && ~error_indicator.identification_moments
         if normalize_jacobians
             norm_dMOMENTS = max(abs(si_dMOMENTS),[],2);
             norm_dMOMENTS = norm_dMOMENTS(:,ones(totparam_nbr,1));
@@ -469,7 +467,7 @@ if info(1) == 0 %no errors in solution
 
         if advanced
             % here we do not normalize (i.e. we set norm_dMOMENTS=1) as the OLS in ident_bruteforce is very sensitive to norm_dMOMENTS
-            [ide_moments.pars, ide_moments.cosndMOMENTS] = ident_bruteforce(dMOMENTS(ind_dMOMENTS,:), max_dim_cova_group, options_.TeX, options_ident.name_tex, options_ident.tittxt, options_ident.tol_deriv);
+            [ide_moments.pars, ide_moments.cosndMOMENTS] = ident_bruteforce(dMOMENTS(ind_dMOMENTS,:), max_dim_cova_group, options_.TeX, options_ident.name_tex, options_ident.tittxt, tol_deriv);
         end
 
         %here we focus on the unnormalized S and V, which is then used in plot_identification.m and for prior_mc > 1
@@ -485,7 +483,7 @@ if info(1) == 0 %no errors in solution
         end
     end
 
-    if ~no_identification_minimal
+    if ~options_ident.no_identification_minimal && ~error_indicator.identification_minimal
         if normalize_jacobians
             ind_dMINIMAL = (find(max(abs(dMINIMAL'),[],1) > tol_deriv)); %index for non-zero rows
             norm_dMINIMAL = max(abs(dMINIMAL(ind_dMINIMAL,:)),[],2);
@@ -499,7 +497,7 @@ if info(1) == 0 %no errors in solution
         ide_minimal.dMINIMAL      = dMINIMAL;
     end
 
-    if ~no_identification_spectrum
+    if ~options_ident.no_identification_spectrum && ~error_indicator.identification_spectrum
         if normalize_jacobians
             ind_dSPECTRUM = (find(max(abs(dSPECTRUM'),[],1) > tol_deriv)); %index for non-zero rows
             tilda_dSPECTRUM = zeros(size(dSPECTRUM));
@@ -523,23 +521,33 @@ if info(1) == 0 %no errors in solution
     if checks_via_subsets
         % identification_checks_via_subsets is only for debugging
         [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = ...
-            identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident);
+            identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident, error_indicator);
+         if ~error_indicator.identification_minimal
+             ide_minimal.minimal_state_space=1;
+         else
+             ide_minimal.minimal_state_space=0;
+         end
     else
         [ide_dynamic.cond, ide_dynamic.rank, ide_dynamic.ind0, ide_dynamic.indno, ide_dynamic.ino, ide_dynamic.Mco, ide_dynamic.Pco, ide_dynamic.jweak, ide_dynamic.jweak_pair] = ...
             identification_checks(dDYNAMIC(ind_dDYNAMIC,:)./norm_dDYNAMIC, 1, tol_rank, tol_sv, modparam_nbr);
-        if ~no_identification_reducedform
+        if ~options_ident.no_identification_reducedform && ~error_indicator.identification_reducedform
             [ide_reducedform.cond, ide_reducedform.rank, ide_reducedform.ind0, ide_reducedform.indno, ide_reducedform.ino, ide_reducedform.Mco, ide_reducedform.Pco, ide_reducedform.jweak, ide_reducedform.jweak_pair] = ...
                 identification_checks(dREDUCEDFORM(ind_dREDUCEDFORM,:)./norm_dREDUCEDFORM, 1, tol_rank, tol_sv, totparam_nbr);
         end
-        if ~no_identification_moments
+        if ~options_ident.no_identification_moments && ~error_indicator.identification_moments
             [ide_moments.cond, ide_moments.rank, ide_moments.ind0, ide_moments.indno, ide_moments.ino, ide_moments.Mco, ide_moments.Pco, ide_moments.jweak, ide_moments.jweak_pair] = ...
                 identification_checks(dMOMENTS(ind_dMOMENTS,:)./norm_dMOMENTS, 1, tol_rank, tol_sv, totparam_nbr);
         end
-        if ~no_identification_minimal
-            [ide_minimal.cond, ide_minimal.rank, ide_minimal.ind0, ide_minimal.indno, ide_minimal.ino, ide_minimal.Mco, ide_minimal.Pco, ide_minimal.jweak, ide_minimal.jweak_pair] = ...
-                identification_checks(dMINIMAL(ind_dMINIMAL,:)./norm_dMINIMAL, 2, tol_rank, tol_sv, totparam_nbr);
+        if ~options_ident.no_identification_minimal 
+            if ~error_indicator.identification_minimal
+                [ide_minimal.cond, ide_minimal.rank, ide_minimal.ind0, ide_minimal.indno, ide_minimal.ino, ide_minimal.Mco, ide_minimal.Pco, ide_minimal.jweak, ide_minimal.jweak_pair] = ...
+                    identification_checks(dMINIMAL(ind_dMINIMAL,:)./norm_dMINIMAL, 2, tol_rank, tol_sv, totparam_nbr);
+                ide_minimal.minimal_state_space=1;
+            else
+                ide_minimal.minimal_state_space=0;
+            end
         end
-        if ~no_identification_spectrum
+        if ~options_ident.no_identification_spectrum && ~error_indicator.identification_spectrum
             [ide_spectrum.cond, ide_spectrum.rank, ide_spectrum.ind0, ide_spectrum.indno, ide_spectrum.ino, ide_spectrum.Mco, ide_spectrum.Pco, ide_spectrum.jweak, ide_spectrum.jweak_pair] = ...
                 identification_checks(tilda_dSPECTRUM, 3, tol_rank, tol_sv, totparam_nbr);
         end
diff --git a/matlab/identification_checks_via_subsets.m b/matlab/identification_checks_via_subsets.m
index d91bee797317865802e06eaeb76841861c4dd38b..1664fded8dc2c2a2f8493d31b07f2e59929ae615 100644
--- a/matlab/identification_checks_via_subsets.m
+++ b/matlab/identification_checks_via_subsets.m
@@ -1,5 +1,5 @@
-function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident)
-%[ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident)
+function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident,error_indicator)
+%[ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident,error_indicator)
 % -------------------------------------------------------------------------
 % Finds problematic sets of paramters via checking the necessary rank condition
 % of the Jacobians for all possible combinations of parameters. The rank is
@@ -23,22 +23,23 @@ function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal]
 % INPUTS
 %   ide_reducedform:    [structure] Containing results from identification
 %                       analysis based on the reduced-form solution (Ratto
-%                       and Iskrev, 2011). If ide_reducedform.no_identification_reducedform
-%                       is 1 then the search for problematic parameter sets will be skipped
+%                       and Iskrev, 2011). If either options_ident.no_identification_reducedform
+%                       or error_indicator.identification_reducedform are 1 then the search for problematic parameter sets will be skipped
 %   ide_moments:        [structure] Containing results from identification
-%                       analysis based on moments (Iskrev, 2010). If
-%                       ide_moments.no_identification_moments is 1 then the search for
+%                       analysis based on moments (Iskrev, 2010). If either
+%                       options_ident.no_identification_moments or error_indicator.identification_moments are 1 then the search for
 %                       problematic parameter sets will be skipped
 %   ide_spectrum:       [structure] Containing results from identification
 %                       analysis based on the spectrum (Qu and Tkachenko, 2012).
-%                       If ide_spectrum.no_identification_spectrum is 1 then the search for
+%                       If either options_ident.no_identification_spectrum or error_indicator.identification_spectrum are 1 then the search for
 %                       problematic parameter sets will be skipped
 %   ide_minimal:        [structure] Containing results from identification
 %                       analysis based on the minimal state space system
-%                       (Komunjer and Ng, 2011). If ide_minimal.no_identification_minimal
-%                       is 1 then the search for problematic parameter sets will be skipped
+%                       (Komunjer and Ng, 2011). If either options_ident.no_identification_minimal or
+%                       error_indicator.identification_minimal are 1 then the search for problematic parameter sets will be skipped
 %   totparam_nbr:       [integer] number of estimated stderr, corr and model parameters
 %   numzerotolrank:     [double] tolerance level for rank compuations
+%   error_indicator     [structure] indicators whether objects could be computed
 % -------------------------------------------------------------------------
 % OUTPUTS
 %   ide_reducedform, ide_moments, ide_spectrum, ide_minimal are augmented by the
@@ -51,7 +52,7 @@ function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal]
 % This function is called by
 %   * identification_analysis.m
 % =========================================================================
-% Copyright (C) 2019 Dynare Team
+% Copyright (C) 2019-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -113,7 +114,7 @@ else
 end
 
 % initialize for reduced form solution criteria
-if ~no_identification_reducedform
+if ~no_identification_reducedform && ~error_indicator.identification_reducedform
     dREDUCEDFORM = ide_reducedform.dREDUCEDFORM;
     dREDUCEDFORM(ide_reducedform.ind_dREDUCEDFORM,:) = dREDUCEDFORM(ide_reducedform.ind_dREDUCEDFORM,:)./ide_reducedform.norm_dREDUCEDFORM; %normalize
     if strcmp(tol_rank,'robust')
@@ -136,7 +137,7 @@ else
 end
 
 % initialize for moments criteria
-if ~no_identification_moments
+if ~no_identification_moments && ~error_indicator.identification_moments
     dMOMENTS = ide_moments.dMOMENTS;
     dMOMENTS(ide_moments.ind_dMOMENTS,:) = dMOMENTS(ide_moments.ind_dMOMENTS,:)./ide_moments.norm_dMOMENTS; %normalize
     if strcmp(tol_rank,'robust')
@@ -159,7 +160,7 @@ else
 end
 
 % initialize for spectrum criteria
-if ~no_identification_spectrum
+if ~no_identification_spectrum && ~error_indicator.identification_spectrum
     dSPECTRUM = ide_spectrum.tilda_dSPECTRUM; %tilda dSPECTRUM is normalized dSPECTRUM matrix in identification_analysis.m
     %alternative normalization
     %dSPECTRUM = ide_spectrum.dSPECTRUM;
@@ -184,7 +185,7 @@ else
 end
 
 % initialize for minimal system criteria
-if ~no_identification_minimal
+if ~no_identification_minimal && ~error_indicator.identification_minimal
     dMINIMAL = ide_minimal.dMINIMAL;
     dMINIMAL(ide_minimal.ind_dMINIMAL,:) = dMINIMAL(ide_minimal.ind_dMINIMAL,:)./ide_minimal.norm_dMINIMAL; %normalize
     dMINIMAL_par  = dMINIMAL(:,1:totparam_nbr);       %part of dMINIMAL that is dependent on parameters
@@ -223,7 +224,7 @@ for j=1:totparam_nbr
             end
         end
     end
-    if ~no_identification_reducedform
+    if ~no_identification_reducedform && ~error_indicator.identification_reducedform
         % Columns correspond to single parameters, i.e. full rank would be equal to 1
         if strcmp(tol_rank,'robust')
             if rank(full(dREDUCEDFORM(:,j))) == 0
@@ -235,7 +236,7 @@ for j=1:totparam_nbr
             end
         end
     end
-    if ~no_identification_moments        
+    if ~no_identification_moments && ~error_indicator.identification_moments
         % Columns correspond to single parameters, i.e. full rank would be equal to 1
         if strcmp(tol_rank,'robust')
             if rank(full(dMOMENTS(:,j))) == 0
@@ -247,13 +248,13 @@ for j=1:totparam_nbr
             end
         end
     end
-    if ~no_identification_spectrum
+    if ~no_identification_spectrum && ~error_indicator.identification_spectrum
         % Diagonal values correspond to single parameters, absolute value needs to be greater than tolerance level
         if abs(dSPECTRUM(j,j)) < tol_rank
             spectrum_problpars{1} = [spectrum_problpars{1};j];
         end
     end
-    if ~no_identification_minimal
+    if ~no_identification_minimal && ~error_indicator.identification_minimal
         % Columns of dMINIMAL_par correspond to single parameters, needs to be augmented with dMINIMAL_rest (part that is independent of parameters),
         % full rank would be equal to 1+dMINIMAL_fixed_rank_nbr
         if strcmp(tol_rank,'robust')
@@ -278,7 +279,7 @@ if ~no_identification_dynamic
         indparam_dDYNAMIC(dynamic_problpars{1}) = []; %remove single unidentified parameters from higher-order sets of indparam
     end
 end
-if ~no_identification_reducedform
+if ~no_identification_reducedform && ~error_indicator.identification_reducedform
     if size(reducedform_problpars{1},1) == (totparam_nbr - rank_dREDUCEDFORM)
         %found all nonidentified parameter sets
         no_identification_reducedform = 1; %skip in the following
@@ -287,7 +288,7 @@ if ~no_identification_reducedform
         indparam_dREDUCEDFORM(reducedform_problpars{1}) = []; %remove single unidentified parameters from higher-order sets of indparam
     end
 end
-if ~no_identification_moments
+if ~no_identification_moments && ~error_indicator.identification_moments
     if size(moments_problpars{1},1) == (totparam_nbr - rank_dMOMENTS)
         %found all nonidentified parameter sets
         no_identification_moments = 1; %skip in the following
@@ -296,7 +297,7 @@ if ~no_identification_moments
         indparam_dMOMENTS(moments_problpars{1}) = []; %remove single unidentified parameters from higher-order sets of indparam
     end
 end
-if ~no_identification_spectrum
+if ~no_identification_spectrum && ~error_indicator.identification_spectrum
     if size(spectrum_problpars{1},1) == (totparam_nbr - rank_dSPECTRUM)
         %found all nonidentified parameter sets
         no_identification_spectrum = 1; %skip in the following
@@ -305,7 +306,7 @@ if ~no_identification_spectrum
         indparam_dSPECTRUM(spectrum_problpars{1}) = []; %remove single unidentified parameters from higher-order sets of indparam
     end
 end
-if ~no_identification_minimal
+if ~no_identification_minimal && ~error_indicator.identification_minimal
     if size(minimal_problpars{1},1) == (totparam_nbr + dMINIMAL_fixed_rank_nbr - rank_dMINIMAL)
         %found all nonidentified parameter sets
         no_identification_minimal = 1; %skip in the following
@@ -323,7 +324,11 @@ indtotparam = unique([indparam_dDYNAMIC indparam_dREDUCEDFORM indparam_dMOMENTS
 for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subsets
     h = dyn_waitbar(0,['Brute force collinearity for ' int2str(j) ' parameters.']);
     %Step1: get all possible unique subsets of j elements
-    if ~no_identification_dynamic || ~no_identification_reducedform || ~no_identification_moments || ~no_identification_spectrum || ~no_identification_minimal
+    if ~no_identification_dynamic ...
+            || (~no_identification_reducedform && ~error_indicator.identification_reducedform)...
+            || (~no_identification_moments && ~error_indicator.identification_moments)...
+            || (~no_identification_spectrum && ~error_indicator.identification_spectrum)...
+            || (~no_identification_minimal && ~error_indicator.identification_minimal)
         indexj=nchoosek(int16(indtotparam),j);  %  int16 speeds up nchoosek
         % One could also use a mex version of nchoosek to speed this up, e.g.VChooseK from https://de.mathworks.com/matlabcentral/fileexchange/26190-vchoosek
     end
@@ -335,25 +340,25 @@ for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subset
     else
         indexj_dDYNAMIC = [];
     end
-    if ~no_identification_reducedform
+    if ~no_identification_reducedform && ~error_indicator.identification_reducedform
         indexj_dREDUCEDFORM = RemoveProblematicParameterSets(indexj,reducedform_problpars);
         rankj_dREDUCEDFORM  = zeros(size(indexj_dREDUCEDFORM,1),1);
     else
         indexj_dREDUCEDFORM = [];
     end
-    if ~no_identification_moments
+    if ~no_identification_moments && ~error_indicator.identification_moments
         indexj_dMOMENTS = RemoveProblematicParameterSets(indexj,moments_problpars);
         rankj_dMOMENTS = zeros(size(indexj_dMOMENTS,1),1);
     else
         indexj_dMOMENTS = [];
     end
-    if ~no_identification_spectrum
+    if ~no_identification_spectrum && ~error_indicator.identification_spectrum
         indexj_dSPECTRUM = RemoveProblematicParameterSets(indexj,spectrum_problpars);
         rankj_dSPECTRUM = zeros(size(indexj_dSPECTRUM,1),1);
     else
         indexj_dSPECTRUM = [];
     end
-    if ~no_identification_minimal
+    if ~no_identification_minimal && ~error_indicator.identification_minimal
         indexj_dMINIMAL = RemoveProblematicParameterSets(indexj,minimal_problpars);
         rankj_dMINIMAL = zeros(size(indexj_dMINIMAL,1),1);
     else
@@ -375,7 +380,7 @@ for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subset
                 end
             end
         end
-        if ~no_identification_reducedform
+        if ~no_identification_reducedform && ~error_indicator.identification_reducedform
             k_dREDUCEDFORM = k_dREDUCEDFORM+1;
             if k_dREDUCEDFORM <= size(indexj_dREDUCEDFORM,1)
                 dREDUCEDFORM_j = dREDUCEDFORM(:,indexj_dREDUCEDFORM(k_dREDUCEDFORM,:)); % pick columns that correspond to parameter subset
@@ -386,7 +391,7 @@ for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subset
                 end
             end
         end
-        if ~no_identification_moments
+        if ~no_identification_moments && ~error_indicator.identification_moments
             k_dMOMENTS = k_dMOMENTS+1;
             if k_dMOMENTS <= size(indexj_dMOMENTS,1)
                 dMOMENTS_j = dMOMENTS(:,indexj_dMOMENTS(k_dMOMENTS,:)); % pick columns that correspond to parameter subset
@@ -397,7 +402,7 @@ for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subset
                 end
             end
         end
-        if ~no_identification_minimal
+        if ~no_identification_minimal && ~error_indicator.identification_minimal
             k_dMINIMAL = k_dMINIMAL+1;
             if k_dMINIMAL <= size(indexj_dMINIMAL,1)
                 dMINIMAL_j = [dMINIMAL_par(:,indexj_dMINIMAL(k_dMINIMAL,:)) dMINIMAL_rest]; % pick columns in dMINIMAL_par that correspond to parameter subset and augment with parameter-indepdendet part dMINIMAL_rest
@@ -408,7 +413,7 @@ for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subset
                 end
             end
         end
-        if ~no_identification_spectrum
+        if ~no_identification_spectrum && ~error_indicator.identification_spectrum
             k_dSPECTRUM = k_dSPECTRUM+1;
             if k_dSPECTRUM <= size(indexj_dSPECTRUM,1)
                 dSPECTRUM_j = dSPECTRUM(indexj_dSPECTRUM(k_dSPECTRUM,:),indexj_dSPECTRUM(k_dSPECTRUM,:)); % pick rows and columns that correspond to parameter subset
@@ -426,16 +431,16 @@ for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subset
     if ~no_identification_dynamic
         dynamic_problpars{j} = indexj_dDYNAMIC(rankj_dDYNAMIC < j,:);
     end
-    if ~no_identification_reducedform
+    if ~no_identification_reducedform && ~error_indicator.identification_reducedform
         reducedform_problpars{j} = indexj_dREDUCEDFORM(rankj_dREDUCEDFORM < j,:);
     end
-    if ~no_identification_moments
+    if ~no_identification_moments && ~error_indicator.identification_moments
         moments_problpars{j} = indexj_dMOMENTS(rankj_dMOMENTS < j,:);
     end
-    if ~no_identification_minimal
+    if ~no_identification_minimal && ~error_indicator.identification_minimal
         minimal_problpars{j} = indexj_dMINIMAL(rankj_dMINIMAL < (j+dMINIMAL_fixed_rank_nbr),:);
     end
-    if ~no_identification_spectrum
+    if ~no_identification_spectrum && ~error_indicator.identification_spectrum
         spectrum_problpars{j} = indexj_dSPECTRUM(rankj_dSPECTRUM < j,:);
     end
 %     % Optional Step 5: % remove redundant 2-sets, eg. if the problematic sets are [(p1,p2);(p1,p3);(p2,p3)], then the unique problematic parameter sets are actually only [(p1,p2),(p1,p3)]