From d43a057af3c7ad301afef3471358e88df2cb5bed Mon Sep 17 00:00:00 2001
From: Marco Ratto <marco.ratto@ec.europa.eu>
Date: Thu, 6 May 2021 16:50:22 +0200
Subject: [PATCH 1/2] store_smoother_results: fix header explanation of output
 arguments

Variances are now in declaration order.
---
 matlab/store_smoother_results.m | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/matlab/store_smoother_results.m b/matlab/store_smoother_results.m
index cc529fab2..645289445 100644
--- a/matlab/store_smoother_results.m
+++ b/matlab/store_smoother_results.m
@@ -34,8 +34,9 @@ function [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,da
 %                   oo_.Smoother.Variance: one-step ahead forecast error variance (declaration order)
 %                   oo_.Smoother.Constant: structure storing the constant term of the smoother
 %                   oo_.Smoother.Trend: structure storing the trend term of the smoother
-%                   oo_.FilteredVariablesKStepAhead: k-step ahead forecast error variance matrices (decision-rule order)
-%                   oo_.FilteredVariablesShockDecomposition: shock decomposition of k-step ahead filtered variables (decision-rule order)
+%                   oo_.FilteredVariablesKStepAhead: k-step ahead filtered variables (declaration order)
+%                   oo_.FilteredVariablesKStepAheadVariances: k-step ahead forecast error variance matrices (declaration order)
+%                   oo_.FilteredVariablesShockDecomposition: shock decomposition of k-step ahead filtered variables (declaration order)
 %                   oo_.FilteredVariables: structure storing the filtered variables
 %                   oo_.UpdatedVariables: structure storing the updated variables
 %                   oo_.SmoothedShocks: structure storing the smoothed shocks

From 551917581fba7858f8f590c508b791716a1ace37 Mon Sep 17 00:00:00 2001
From: Marco Ratto <marco.ratto@ec.europa.eu>
Date: Thu, 6 May 2021 16:56:03 +0200
Subject: [PATCH 2/2] Bug fixes in fitered variances of smoother

- kalman_algo=1: kstep-ahead variances were WRONG, since Pf was initialized using P in PREVIOUS period
- kalman_algo=2: output argument for filtered varainces should be P1, not P (P are UPDATED variances, there).

For kalman_algo=2, also make a small factorization fix (compute P(:,:,t+1) before defining Pf, so to compute 1-step ahead variance only once)
---
 matlab/missing_DiffuseKalmanSmootherH1_Z.m |  6 ++++--
 matlab/missing_DiffuseKalmanSmootherH3_Z.m | 15 +++++++++------
 2 files changed, 13 insertions(+), 8 deletions(-)

diff --git a/matlab/missing_DiffuseKalmanSmootherH1_Z.m b/matlab/missing_DiffuseKalmanSmootherH1_Z.m
index 2a4ed39c4..fbdf7e211 100644
--- a/matlab/missing_DiffuseKalmanSmootherH1_Z.m
+++ b/matlab/missing_DiffuseKalmanSmootherH1_Z.m
@@ -211,10 +211,12 @@ while notsteady && t<smpl
         P(:,:,t+1)  = T*P(:,:,t)*L(:,:,t)' + QQ;
     end
     a(:,t+1)    = T*atilde(:,t);
-    Pf          = P(:,:,t);
+    Pf          = P(:,:,t+1);
     aK(1,:,t+1) = a(:,t+1);
     for jnk=1:nk
-        Pf = T*Pf*T' + QQ;
+        if jnk>1
+            Pf = T*Pf*T' + QQ;
+        end
         PK(jnk,:,:,t+jnk) = Pf;
         if jnk>1
             aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
diff --git a/matlab/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
index cfdbdd08c..ff30820ef 100644
--- a/matlab/missing_DiffuseKalmanSmootherH3_Z.m
+++ b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
@@ -1,5 +1,5 @@
-function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp,V] = missing_DiffuseKalmanSmootherH3_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag)
-% function [alphahat,epsilonhat,etahat,a1,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag,state_uncertainty_flag)
+function [alphahat,epsilonhat,etahat,a,P1,aK,PK,decomp,V] = missing_DiffuseKalmanSmootherH3_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag)
+% function [alphahat,epsilonhat,etahat,a1,P1,aK,PK,d,decomp] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag,state_uncertainty_flag)
 % Computes the diffuse kalman smoother in the case of a singular var-cov matrix.
 % Univariate treatment of multivariate time series.
 %
@@ -31,7 +31,7 @@ function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp,V] = missing_DiffuseKalman
 %    a:        matrix of updated variables (a_{t|t})
 %    aK:       3D array of k step ahead filtered state variables (a_{t+k|t})
 %              (meaningless for periods 1:d)
-%    P:        3D array of one-step ahead forecast error variance
+%    P1:        3D array of one-step ahead forecast error variance
 %              matrices
 %    PK:       4D array of k-step ahead forecast error variance
 %              matrices (meaningless for periods 1:d)
@@ -221,18 +221,21 @@ while notsteady && t<smpl
         end
     end
     a1(:,t+1) = T*a(:,t);                                                   %transition according to (6.14) in DK (2012)
-    Pf          = P(:,:,t);
+    P(:,:,t+1) = T*P(:,:,t)*T' + QQ;                                        %transition according to (6.14) in DK (2012)
+    Pf          = P(:,:,t+1);
     aK(1,:,t+1) = a1(:,t+1);
     for jnk=1:nk
-        Pf = T*Pf*T' + QQ;
+        if jnk>1
+            Pf = T*Pf*T' + QQ;
+        end
         PK(jnk,:,:,t+jnk) = Pf;
         if jnk>1
             aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
         end
     end
-    P(:,:,t+1) = T*P(:,:,t)*T' + QQ;                                        %transition according to (6.14) in DK (2012)
                                                                             %  notsteady   = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<kalman_tol);
 end
+P1(:,:,t+1) = P(:,:,t+1);
 % $$$ P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
 % $$$ P1_s=tril(P1(:,:,t))+tril(P1(:,:,t),-1)';
 % $$$ Fi_s = Fi(:,t);