diff --git a/matlab/+occbin/solve_one_constraint.m b/matlab/+occbin/solve_one_constraint.m
index 27d01a720..5b710cff5 100644
--- a/matlab/+occbin/solve_one_constraint.m
+++ b/matlab/+occbin/solve_one_constraint.m
@@ -125,20 +125,24 @@ for shock_period = 1:n_shocks_periods
     
     regime_change_this_iteration=true;
     iter = 0;
+    nperiods_endogenously_increased = false;
     guess_history_it = false;
     if guess_history && (shock_period<=length(regime_history_guess)) %beyond guess regime history
         guess_history_it = true;
     end
     
     is_periodic=false;
+    is_periodic_loop=false;
     binding_indicator_history={};
     max_err = NaN(max_iter,1);
+    regime_violates_constraint_in_expectation = false(max_iter,1);
     
-    while (regime_change_this_iteration && iter<max_iter && ~is_periodic)
+    while (regime_change_this_iteration && iter<max_iter && ~is_periodic && ~is_periodic_loop)
         iter = iter +1;
         if binding_indicator(end) && nperiods_0<opts_simul_.max_periods
             binding_indicator = [binding_indicator; false ];
             nperiods_0 = nperiods_0 + 1;
+            nperiods_endogenously_increased = true;
             disp_verbose(['nperiods has been endogenously increased up to ' int2str(nperiods_0) '.'],opts_simul_.debug)
         elseif nperiods_0>=opts_simul_.max_periods
             % enforce endogenously increased nperiods to not violate max_periods
@@ -201,6 +205,7 @@ for shock_period = 1:n_shocks_periods
                 err_relax = err.relax_constraint_1(relax.constraint_1(1:end_periods) & binding_indicator(1:end_periods));
                 max_err(iter) = max(abs([err_viol;err_relax]));
                 regime_change_this_iteration = true;
+                regime_violates_constraint_in_expectation(iter) = any(binding.constraint_1(1:end_periods) & ~binding_indicator(1:end_periods));
             else
                 regime_change_this_iteration = false;
                 max_err(iter) = 0;
@@ -220,16 +225,43 @@ for shock_period = 1:n_shocks_periods
                 binding_indicator= (binding_indicator | binding.constraint_1) & ~(binding_indicator & relax.constraint_1);
             end
             
-            if iter>1 && regime_change_this_iteration
+            if iter>1 && regime_change_this_iteration && ~nperiods_endogenously_increased
+                % check for periodic solution only if nperiods is not
+                % increased endogenously
+                % first check for infinite loop
+                is_periodic_loop = false(iter-1,1);
                 for kiter=1:iter-1
-                    vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)];
-                    is_periodic(kiter) = isequal(vvv, binding_indicator);
+                    if size(binding_indicator,1)== size(binding_indicator_history{kiter},1)
+                        %                     vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)];
+                        %                     is_periodic(kiter) = isequal(vvv, binding_indicator);
+                        is_periodic_loop(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator);
+                    else
+                        is_periodic_loop(kiter) = false;
+                    end
+                end
+                is_periodic_loop_all =is_periodic_loop;
+                is_periodic_loop =  any(is_periodic_loop);
+                % only accept periodicity where regimes differ by one
+                % period!
+                is_periodic = false(iter-1,1);
+                for kiter=iter-1
+                    if size(binding_indicator,1)== size(binding_indicator_history{kiter},1)
+                        %                     vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)];
+                        %                     is_periodic(kiter) = isequal(vvv, binding_indicator);
+                        is_periodic(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator) && length(find(binding_indicator_history{iter}-binding_indicator))==1;
+                    else
+                        is_periodic(kiter)=false;
+                    end
                 end
                 is_periodic_all =is_periodic;
                 is_periodic =  any(is_periodic);
                 if is_periodic && periodic_solution
-                    [merr,imerr]=min(max_err(find(is_periodic_all,1):end));
                     inx = find(is_periodic_all,1):iter;
+                    inx1 = inx(find(~regime_violates_constraint_in_expectation(inx)));
+                    if not(isempty(inx1))
+                        inx=inx1;
+                    end
+                    [merr,imerr]=min(max_err(inx));
                     inx = inx(imerr);
                     binding_indicator=binding_indicator_history{inx};
                     if inx<iter
@@ -281,9 +313,14 @@ for shock_period = 1:n_shocks_periods
                     return
                 end
             else
+                if is_periodic_loop
+                    disp_verbose('Did not converge -- infinite loop of regimes.',opts_simul_.debug)
+                    error_flag = 313;
+                else
                 disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
-                if opts_simul_.waitbar; dyn_waitbar_close(hh); end
                     error_flag = 311;
+                end
+                if opts_simul_.waitbar; dyn_waitbar_close(hh); end
                 return
             end
         else
diff --git a/matlab/+occbin/solve_two_constraints.m b/matlab/+occbin/solve_two_constraints.m
index 4e448b241..15c8ae1c2 100644
--- a/matlab/+occbin/solve_two_constraints.m
+++ b/matlab/+occbin/solve_two_constraints.m
@@ -133,20 +133,24 @@ for shock_period = 1:n_shocks_periods
         dyn_waitbar(shock_period/n_shocks_periods, hh, sprintf('Period %u of %u', shock_period,n_shocks_periods));
     end
     regime_change_this_iteration=true;
+    nperiods_endogenously_increased = false;
     iter = 0;
     guess_history_it = false;
     if guess_history && (shock_period<=length(regime_history_guess)) %beyond guess regime history
         guess_history_it = true;
     end
     is_periodic=false;
+    is_periodic_loop=false;
     binding_indicator_history={};
     max_err = NaN(max_iter,1);
+    regime_violates_constraint_in_expectation = false(max_iter,1);
     
-    while (regime_change_this_iteration && iter<max_iter && ~is_periodic)
+    while (regime_change_this_iteration && iter<max_iter && ~is_periodic && ~is_periodic_loop)
         iter = iter +1;
         if any(binding_indicator(end,:)) && nperiods_0<opts_simul_.max_periods
             binding_indicator = [binding_indicator; false(1,2)];
             nperiods_0 = nperiods_0 + 1;
+            nperiods_endogenously_increased = true;
             disp_verbose(['nperiods has been endogenously increased up to ' int2str(nperiods_0) '.'],opts_simul_.debug)
         elseif nperiods_0>=opts_simul_.max_periods
             % enforce endogenously increased nperiods to not violate max_periods
@@ -248,17 +252,43 @@ for shock_period = 1:n_shocks_periods
             end
             binding_indicator = reshape(binding_indicator,nperiods_0+1,2);
             
-            if iter>1 && regime_change_this_iteration
-                is_periodic=false(1,iter-1);
+            if iter>1 && regime_change_this_iteration && ~nperiods_endogenously_increased
+                % check for periodic solution only if nperiods is not
+                % increased endogenously
+                % first check for infinite loop
+                is_periodic_loop = false(iter-1,1);
                 for kiter=1:iter-1
-                    vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 2)];
-                    is_periodic(kiter) = isequal(vvv, binding_indicator);
+                    if size(binding_indicator,1)== size(binding_indicator_history{kiter},1)
+                        %                     vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)];
+                        %                     is_periodic(kiter) = isequal(vvv, binding_indicator);
+                        is_periodic_loop(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator);
+                    else
+                        is_periodic_loop(kiter) = false;
+                    end
+                end
+%                 is_periodic_loop_all =is_periodic_loop;
+                is_periodic_loop =  any(is_periodic_loop);
+                % only accept periodicity where regimes differ by one
+                % period!
+                is_periodic=false(1,iter-1);
+                for kiter=iter-1
+                    if size(binding_indicator,1)== size(binding_indicator_history{kiter},1)
+                        %                     vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)];
+                        %                     is_periodic(kiter) = isequal(vvv, binding_indicator);
+                        is_periodic(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator) && length(find(binding_indicator_history{iter}(:,1)-binding_indicator(:,1)))<=1  && length(find(binding_indicator_history{iter}(:,2)-binding_indicator(:,2)))<=1;
+                    else
+                        is_periodic(kiter)=false;
+                    end
                 end
                 is_periodic_all = is_periodic;
                 is_periodic =  any(is_periodic);
                 if is_periodic && periodic_solution
-                    [min_err,index_min_err]=min(max_err(find(is_periodic_all,1):end));
                     inx = find(is_periodic_all,1):iter;
+                    inx1 = inx(find(~regime_violates_constraint_in_expectation(inx)));
+                    if not(isempty(inx1))
+                        inx=inx1;
+                    end
+                    [min_err,index_min_err]=min(max_err(inx));
                     inx = inx(index_min_err);
                     binding_indicator=binding_indicator_history{inx}; %select regime history with same result, but smallest error
                     if inx<iter %update if needed
@@ -312,8 +342,13 @@ for shock_period = 1:n_shocks_periods
                     return;
                 end
             else
+                if is_periodic_loop
+                    disp_verbose('Did not converge -- infinite loop of regimes.',opts_simul_.debug)
+                    error_flag = 313;
+                else
                 disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
                 error_flag = 311;
+                end
                 if opts_simul_.waitbar; dyn_waitbar_close(hh); end
                 return;
             end
diff --git a/matlab/get_error_message.m b/matlab/get_error_message.m
index 8b463f619..9d9cdae45 100644
--- a/matlab/get_error_message.m
+++ b/matlab/get_error_message.m
@@ -188,6 +188,8 @@ switch info(1)
         message = 'Occbin: Simulation did not converge, increase maxit or check_ahead_periods.';        
     case 312
         message = 'Occbin: Constraint(s) are binding at the end of the sample.';        
+    case 313
+        message = 'Occbin: Simulation did not converge -- infinite loop of regimes';        
     case 320
         message = 'Piecewise linear Kalman filter: There was a problem in obtaining the likelihood.';
     case 401