diff --git a/matlab/DsgeLikelihood.m b/matlab/DsgeLikelihood.m
index fb0ff5177..e30068fbd 100644
--- a/matlab/DsgeLikelihood.m
+++ b/matlab/DsgeLikelihood.m
@@ -170,15 +170,7 @@ if options_.lik_init == 1               % Kalman filter
     if kalman_algo ~= 2
         kalman_algo = 1;
     end
-    [Pstar,junk,unit_roots] = lyapunov_symm(T,R*Q*R',options_.qz_criterium,...
-                                            options_.lyapunov_complex_threshold);
-    if ~isempty(unit_roots)
-        % if unit roots the penalty equals the sum of distance to 2-qz_criterium
-        fval = bayestopt_.penalty + sum(unit_roots-2+ ...
-                                        options_.qz_criterium);
-        cost_flag = 0;
-        return
-    end
+    Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);
     Pinf        = [];
 elseif options_.lik_init == 2   % Old Diffuse Kalman filter
     if kalman_algo ~= 2
diff --git a/matlab/lyapunov_symm.m b/matlab/lyapunov_symm.m
index 5cf180daf..801df633f 100644
--- a/matlab/lyapunov_symm.m
+++ b/matlab/lyapunov_symm.m
@@ -1,4 +1,4 @@
-function [x,u,unit_roots] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,method)
+function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,method)
 % Solves the Lyapunov equation x-a*x*a' = b, for b and x symmetric matrices.
 % If a has some unit roots, the function computes only the solution of the stable subsystem.
 %  
@@ -15,7 +15,6 @@ function [x,u,unit_roots] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_thre
 % OUTPUTS
 %   x:      [double]    m*m solution matrix of the lyapunov equation, where m is the dimension of the stable subsystem.
 %   u:      [double]    Schur vectors associated with unit roots  
-%   unit_roots [double] vector containing roots too close to 1
 %
 % ALGORITHM
 %   Uses reordered Schur decomposition
@@ -58,13 +57,10 @@ if size(a,1) == 1
     return
 end
 
-unit_roots = [];
 if method<2
     [U,T] = schur(a);
-    roots = abs(ordeig(T));
-    e1 = roots > 2-qz_criterium;
-    k = sum(e1);       % Number of unit roots.
-    unit_roots = roots(1:k);
+    e1 = abs(ordeig(T)) > 2-qz_criterium;
+    k = sum(e1);       % Number of unit roots. 
     n = length(e1)-k;  % Number of stationary variables.
     if k > 0
         % Selects stable roots