diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index ee92a99031b3f4872451bc04766fa3f0a6e09d71..abc46dbe8cf26a9cb60e1f409efd171ea5758ffc 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -319,7 +319,7 @@ if kalman_algo == 2 || kalman_algo == 4
     a_initial     = zeros(np,1);
     a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_);
     a_initial=ST*a_initial; %set state prediction for first Kalman step;
-    [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, regimes_,TT,RR,CC] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ...
+    [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, regimes_,TT,RR,CC,TTx,RRx,CCx] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ...
         Z,R1,Q,diag(H), ...
         Pinf,Pstar,data1,vobs,np,smpl,data_index, ...
         options_.nk,kalman_tol,diffuse_kalman_tol, ...
@@ -386,9 +386,11 @@ else
         else
             opts_simul = options_.occbin.simul;
         end
+        % reconstruct smoothed variables
         aaa=zeros(M_.endo_nbr,gend);
         aaa(oo_.dr.restrict_var_list,:)=alphahat;
-        for k=2:gend
+        iTx = zeros(size(TTx));
+        for k=1:gend
             if isoccbin
                 A = TT(:,:,k);
                 B = RR(:,:,k);
@@ -396,20 +398,53 @@ else
             else
                 C=0;
             end
-            aaa(:,k) = C+A*aaa(:,k-1)+B*etahat(:,k);
+            iT = pinv(TTx(:,:,k));
+            % store pinv
+            iTx(:,:,k) = iT;
+            Tstar = A(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),oo_.dr.restrict_var_list);
+            Rstar = B(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),:);
+            Cstar = C(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list));
+            AS = Tstar*iT;
+            BS = Rstar-AS*RRx(:,:,k);
+            CS = Cstar-AS*CCx(:,k);
+            static_var_list = ~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list);
+            ilagged = any(abs(AS*TTx(:,:,k)-Tstar)'>1.e-12);
+            static_var_list0 = static_var_list;
+            static_var_list0(static_var_list) = ilagged;
+            static_var_list(static_var_list) = ~ilagged;
+            aaa(static_var_list,k) = AS(~ilagged,:)*alphahat(:,k)+BS(~ilagged,:)*etahat(:,k)+CS(~ilagged);
+            if any(ilagged) && k>1
+                aaa(static_var_list0,k) = Tstar(ilagged,:)*alphahat(:,k-1)+Rstar(ilagged,:)*etahat(:,k)+Cstar(ilagged);
+            end
+            
         end
         alphahat=aaa;
-        aaa=zeros(M_.endo_nbr,gend);
+
+        % reconstruct updated variables
         bbb=zeros(M_.endo_nbr,gend);
-        bbb(oo_.dr.restrict_var_list,:)=ahat;
-        aaa(oo_.dr.restrict_var_list,:)=aahat;
-        for k=d+2:gend
+        bbb(oo_.dr.restrict_var_list,:)=ahat; % this is t|t
+        for k=1:gend
             if isoccbin
                 A = TT(:,:,k);
                 B = RR(:,:,k);
                 C = CC(:,k);
-                bbb(:,k) = C+A*aaa(:,k-1)+B*eehat(:,k);
-            else
+                iT = iTx(:,:,k);
+                Tstar = A(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),oo_.dr.restrict_var_list);
+                Rstar = B(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),:);
+                Cstar = C(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list));
+                AS = Tstar*iT;
+                BS = Rstar-AS*RRx(:,:,k);
+                CS = Cstar-AS*CCx(:,k);
+                static_var_list = ~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list);
+                ilagged = any(abs(AS*TTx(:,:,k)-Tstar)'>1.e-12);
+                static_var_list0 = static_var_list;
+                static_var_list0(static_var_list) = ilagged;
+                static_var_list(static_var_list) = ~ilagged;
+                bbb(static_var_list,k) = AS(~ilagged,:)*ahat(:,k)+BS(~ilagged,:)*eehat(:,k)+CS(~ilagged);
+                if any(ilagged) && k>d+1
+                    bbb(static_var_list0,k) = Tstar(ilagged,:)*aahat(:,k-1)+Rstar(ilagged,:)*eehat(:,k)+Cstar(ilagged);
+                end
+            elseif k>d+1
                 opts_simul.curb_retrench = options_.occbin.smoother.curb_retrench;
                 opts_simul.waitbar = options_.occbin.smoother.waitbar;
                 opts_simul.maxit = options_.occbin.smoother.maxit;
@@ -433,8 +468,6 @@ else
                 bbb(oo_.dr.inv_order_var,k) = out.piecewise(1,:) - out.ys';
             end
         end
-        % do not overwrite accurate computations using reduced st. space
-        bbb(oo_.dr.restrict_var_list,:)=ahat;
         ahat0=ahat;
         ahat=bbb;
         if ~isempty(P)
diff --git a/matlab/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
index d585d3ca4eae6045694abd4ba29250c72de5d721..ab9eb3146bae730de8a92e54b2e52dd6b0aad939 100644
--- a/matlab/missing_DiffuseKalmanSmootherH3_Z.m
+++ b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
@@ -486,6 +486,9 @@ varargout{1} = regimes_;
 varargout{2} = TTT;
 varargout{3} = RRR;
 varargout{4} = CCC;
+varargout{5} = TT;
+varargout{6} = RR;
+varargout{7} = CC;
 % $$$ P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
 % $$$ P1_s=tril(P1(:,:,t))+tril(P1(:,:,t),-1)';
 % $$$ Fi_s = Fi(:,t);