diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 6466ad2f1ef6fea5023645365987432abd775b6b..d83aa07d7bbb4571cc5e16b678af722c1c6286cd 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -183,7 +183,7 @@ elseif options_.lik_init == 3           % Diffuse Kalman filter
             Z   = [Z, eye(vobs)];
         end
     end
-    [Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,options_.qz_criterium,oo_.dr.restrict_var_list);
+    [Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,options_.qz_criterium);
 elseif options_.lik_init == 4           % Start from the solution of the Riccati equation.
     [err, Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,vobs)),H);
     mexErrCheck('kalman_steady_state',err);
diff --git a/matlab/compute_Pinf_Pstar.m b/matlab/compute_Pinf_Pstar.m
index d847ffacc28b98ca806c817401085bedd758502b..4eb5fbc87a8bfa7ec5161079e236b0e38d2e2de1 100644
--- a/matlab/compute_Pinf_Pstar.m
+++ b/matlab/compute_Pinf_Pstar.m
@@ -48,8 +48,24 @@ function [Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,qz_criterium, restrict_colum
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 np = size(T,1);
+if nargin == 6
+    indx = restrict_columns;
+    indx0=find(~ismember([1:np],indx));
+else
+%     indx=(find(max(abs(T))>0));
+%     indx0=(find(max(abs(T))==0));
+    indx=(find(max(abs(T))>=1.e-10));
+    indx0=(find(max(abs(T))<1.e-10));
+end
+np0=length(indx0);
+Tbkp = T;
+T0=T(indx0,indx); %static variables vs. dynamic ones
+R0=R(indx0,:); % matrix of shocks for static variables
 
-% perform Kitagawa transformation
+% perform Kitagawa transformation only for non-zero columns of T
+T=T(indx,indx);
+R=R(indx,:);
+np = size(T,1);
 [QT,ST] = schur(T);
 e1 = abs(ordeig(ST)) > 2-qz_criterium;
 [QT,ST] = ordschur(QT,ST,e1);
@@ -59,7 +75,6 @@ nk1 = nk+1;
 Pstar = zeros(np,np);
 R1 = QT'*R;
 B = R1*Q*R1';
-% computes variance of stationary block (lower right)
 i = np;
 while i >= nk+2
     if ST(i,i-1) == 0
@@ -100,13 +115,71 @@ if i == nk+1
     Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1));
 end
 
+if np0
+    ST1=ST;
+    
+    % now I recover stationarized static variables
+    % using 
+    % ss = s-A*z and
+    % z-z(-1) (growth rates of unit roots) only depends on stationary variables
+    Pstar = blkdiag(zeros(np0),Pstar);
+    ST = [zeros(length(Pstar),length(indx0)) [T0*QT ;ST]];
+    R1 = [R0; R1];
+    % build the matrix for stationarized variables
+    STinf = ST(np0+1:np0+nk,np0+1:np0+nk);
+    iSTinf = inv(STinf);
+    ST0=ST;
+    ST0(:,1:np0+nk)=0;  % stationarized static + 1st difference only respond to lagged stationary states
+    ST00 = ST(1:np0,np0+1:np0+nk);
+    %     A\B is the matrix division of A into B, which is roughly the
+    %     same as INV(A)*B 
+    ST0(1:np0,np0+nk+1:end) = ST(1:np0,np0+nk+1:end)-ST00*(iSTinf*ST(np0+1:np0+nk,np0+nk+1:end)); % snip non-stationary part
+    R10 = R1;
+    R10(1:np0,:) = R1(1:np0,:)-ST00*(iSTinf*R1(np0+1:np0+nk,:)); % snip non-stationary part
+    
+    % kill non-stationary part before projecting Pstar
+    ST0(np0+1:np0+nk,:)=0; 
+    R10(np0+1:np0+nk,:)=0; % is this questionable???? IT HAS TO in order to match Michel's version!!!
+
+    % project Pstar onto static x
+    Pstar = ...
+        ST0*Pstar*ST0'+ ...
+        R10*Q*R10';
+    
+    % QT(1:np0,np0+1:np0+nk) = QT(1:np0,np0+1:np0+nk)+ST(1:np0,np0+1:np0+nk);  %%% is this questionable ????
+    % reorder QT entries
+else
+    STinf = ST(np0+1:np0+nk,np0+1:np0+nk);
+
+end
+
 % stochastic trends with no influence on observed variables are
 % arbitrarily initialized to zero
 Pinf = zeros(np,np);
 Pinf(1:nk,1:nk) = eye(nk);
-for k = 1:nk
-    if norm(QT(mf,:)*ST(:,k)) < 1e-8
-        Pinf(k,k) = 0;
+if np0
+    STtriu = STinf-eye(nk);
+%     A\B is the matrix division of A into B, which is roughly the
+%     same as INV(A)*B 
+    STinf0 = ST00*(eye(nk)-iSTinf*STtriu);
+    Pinf = blkdiag(zeros(np0),Pinf);
+    QT = blkdiag(eye(np0),QT);
+    QTinf = QT;
+    QTinf(1:np0,np0+1:np0+nk) = STinf0;
+    QTinf([indx0(:); indx(:)],:) = QTinf;
+    STinf1 = [zeros(np0+np,np0) [STinf0; eye(nk); zeros(np-nk,nk)] zeros(np0+np,np-nk)];
+    for k = 1:nk
+        if norm(QTinf(mf,:)*ST([indx0(:); indx(:)],k+np0)) < 1e-8
+            Pinf(k+np0,k+np0) = 0;
+        end
+    end
+    Pinf = STinf1*Pinf*STinf1';
+    QT([indx0(:); indx(:)],:) = QT;
+else
+    for k = 1:nk
+        if norm(QT(mf,:)*ST(:,k)) < 1e-8
+            Pinf(k+np0,k+np0) = 0;
+        end
     end
 end
 
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index c3e69b0ef1ee80663e4b0f190ae0fbcb85da486f..a6b9071fcb81fa31f7633e2f1e04dfca2693a9a5 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -372,7 +372,7 @@ switch DynareOptions.lik_init
         error(['The model requires Diffuse filter, but you specified a different Kalman filter. You must set options_.kalman_algo ' ...
                'to 0 (default), 3 or 4'])
     end
-    [Pstar,Pinf] = compute_Pinf_Pstar(Z,T,R,Q,DynareOptions.qz_criterium,[1:length(T)]);
+    [Pstar,Pinf] = compute_Pinf_Pstar(Z,T,R,Q,DynareOptions.qz_criterium);
     Z =zeros(length(BayesInfo.mf),size(T,1));
     for i = 1:length(BayesInfo.mf)
         Z(i,BayesInfo.mf(i))=1;