From 570d7de7771d9a1787a27295b114c060d598a43e Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Fri, 26 Jun 2020 18:24:33 +0200
Subject: [PATCH] Fix moment computation with Measurement errors

- check logic for M_.H was faulty
- M_.H was not updated in posterior sampling

(manually cherry picked from commit 6e06acc7f4fd8fed9f8dfc647a021770e952ed15)
---
 matlab/compute_moments_varendo.m                |  4 ++--
 matlab/conditional_variance_decomposition.m     |  2 +-
 matlab/disp_th_moments.m                        |  2 +-
 ...retical_conditional_variance_decomposition.m | 17 +++++++++--------
 matlab/dsge_simulated_theoretical_correlation.m | 14 +++++++-------
 matlab/dsge_simulated_theoretical_covariance.m  | 14 +++++++-------
 ...mulated_theoretical_variance_decomposition.m | 17 +++++++++--------
 matlab/posterior_analysis.m                     |  8 ++++----
 matlab/prior_analysis.m                         |  2 +-
 matlab/prior_sampler.m                          |  2 +-
 matlab/selec_posterior_draws.m                  |  6 +++---
 matlab/stoch_simul.m                            |  2 +-
 12 files changed, 46 insertions(+), 44 deletions(-)

diff --git a/matlab/compute_moments_varendo.m b/matlab/compute_moments_varendo.m
index 4bfabce340..75d6a69da0 100644
--- a/matlab/compute_moments_varendo.m
+++ b/matlab/compute_moments_varendo.m
@@ -144,7 +144,7 @@ if M_.exo_nbr > 1
         skipline();
     end
     skipline();
-    if ~all(M_.H==0)
+    if ~all(diag(M_.H)==0)
         if isoctave || matlab_ver_less_than('8.1')
             [observable_name_requested_vars, varlist_pos] = intersect_stable(var_list_, options_.varobs);
         else
@@ -231,7 +231,7 @@ if M_.exo_nbr > 1
             end
         end
         skipline();
-        if ~all(M_.H==0)
+        if ~all(diag(M_.H)==0)
             if ~isempty(observable_name_requested_vars)
                 NumberOfObservedEndogenousVariables = length(observable_name_requested_vars);
                 temp=NaN(NumberOfObservedEndogenousVariables,NumberOfExogenousVariables+1,length(Steps));
diff --git a/matlab/conditional_variance_decomposition.m b/matlab/conditional_variance_decomposition.m
index 4ef81da553..a77d3ea12a 100644
--- a/matlab/conditional_variance_decomposition.m
+++ b/matlab/conditional_variance_decomposition.m
@@ -88,7 +88,7 @@ end
 % get intersection of requested variables and observed variables with
 % Measurement error
 
-if ~all(StateSpaceModel.measurement_error==0)
+if ~all(diag(StateSpaceModel.measurement_error)==0)
     if isoctave || matlab_ver_less_than('8.1')
         [observable_pos,index_subset,index_observables]=intersect_stable(SubsetOfVariables,StateSpaceModel.observable_pos);
     else
diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m
index 33fdb7100b..3bcaa538eb 100644
--- a/matlab/disp_th_moments.m
+++ b/matlab/disp_th_moments.m
@@ -54,7 +54,7 @@ oo_.mean = m;
 oo_.var = oo_.gamma_y{1};
 
 ME_present=0;
-if ~all(M_.H==0)
+if ~all(diag(M_.H)==0)
     if isoctave || matlab_ver_less_than('8.1')
         [observable_pos_requested_vars,index_subset,index_observables]=intersect_stable(ivar,options_.varobs_id);
     else
diff --git a/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m b/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m
index 540b2e2392..3e20717df9 100644
--- a/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m
+++ b/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m
@@ -82,7 +82,7 @@ NumberOfSavedElementsPerSimulation = nvar*M_.exo_nbr*length(Steps);
 MaXNumberOfConditionalDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation/8);
 
 ME_present=0;
-if ~all(M_.H==0)
+if ~all(diag(M_.H)==0)
     if isoctave || matlab_ver_less_than('8.1')
         [observable_pos_requested_vars,index_subset,index_observables]=intersect_stable(ivar,options_.varobs_id);
     else
@@ -131,20 +131,20 @@ linea = 0;
 linea_ME = 0;
 for file = 1:NumberOfDrawsFiles
     if posterior
-        load([M_.dname '/metropolis/' M_.fname '_' type '_draws' num2str(file) ]);
+        temp=load([M_.dname '/metropolis/' M_.fname '_' type '_draws' num2str(file) ]);
     else
-        load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
+        temp=load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
     end
-    isdrsaved = columns(pdraws)-1;
-    NumberOfDraws = rows(pdraws);
+    isdrsaved = columns(temp.pdraws)-1;
+    NumberOfDraws = rows(temp.pdraws);
     for linee = 1:NumberOfDraws
         linea = linea+1;
         linea_ME = linea_ME+1;
         if isdrsaved
-            M_=set_parameters_locally(M_,pdraws{linee,1});% Needed to update the covariance matrix of the state innovations.
-            dr = pdraws{linee,2};
+            M_=set_parameters_locally(M_,temp.pdraws{linee,1});% Needed to update the covariance matrix of the state innovations.
+            dr = temp.pdraws{linee,2};
         else
-            M_=set_parameters_locally(M_,pdraws{linee,1});
+            M_=set_parameters_locally(M_,temp.pdraws{linee,1});
             [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
         end
         if first_call
@@ -163,6 +163,7 @@ for file = 1:NumberOfDrawsFiles
         end
         [StateSpaceModel.transition_matrix,StateSpaceModel.impulse_matrix] = kalman_transition_matrix(dr,iv,ic,M_.exo_nbr);
         StateSpaceModel.state_innovations_covariance_matrix = M_.Sigma_e;
+        M_ = set_measurement_errors(temp.pdraws{linee,1},temp.estim_params_,M_);
         StateSpaceModel.measurement_error=M_.H;
         clear('dr');
         [ConditionalVarianceDecomposition, ConditionalVarianceDecomposition_ME]=conditional_variance_decomposition(StateSpaceModel, Steps, ivar);
diff --git a/matlab/dsge_simulated_theoretical_correlation.m b/matlab/dsge_simulated_theoretical_correlation.m
index fb00d5ddcd..c32aa9db42 100644
--- a/matlab/dsge_simulated_theoretical_correlation.m
+++ b/matlab/dsge_simulated_theoretical_correlation.m
@@ -94,19 +94,19 @@ CorrFileNumber = 1;
 linea = 0;
 for file = 1:NumberOfDrawsFiles
     if posterior
-        load([M_.dname '/metropolis/' M_.fname '_' type '_draws' num2str(file) ]);
+        temp=load([M_.dname '/metropolis/' M_.fname '_' type '_draws' num2str(file) ]);
     else
-        load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
+        temp=load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
     end
-    NumberOfDraws = rows(pdraws);
-    isdrsaved = columns(pdraws)-1;
+    NumberOfDraws = rows(temp.pdraws);
+    isdrsaved = columns(temp.pdraws)-1;
     for linee = 1:NumberOfDraws
         linea = linea+1;
         if isdrsaved
-            M_=set_parameters_locally(M_,pdraws{linee,1});% Needed to update the covariance matrix of the state innovations.
-            dr = pdraws{linee,2};
+            M_=set_parameters_locally(M_,temp.pdraws{linee,1});% Needed to update the covariance matrix of the state innovations.
+            dr = temp.pdraws{linee,2};
         else
-            M_=set_parameters_locally(M_,pdraws{linee,1});
+            M_=set_parameters_locally(M_,temp.pdraws{linee,1});
             [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
         end
         tmp = th_autocovariances(dr,ivar,M_,options_,nodecomposition);
diff --git a/matlab/dsge_simulated_theoretical_covariance.m b/matlab/dsge_simulated_theoretical_covariance.m
index 0be948bc3c..bbd4702751 100644
--- a/matlab/dsge_simulated_theoretical_covariance.m
+++ b/matlab/dsge_simulated_theoretical_covariance.m
@@ -93,19 +93,19 @@ CovarFileNumber = 1;
 linea = 0;
 for file = 1:NumberOfDrawsFiles
     if posterior
-        load([M_.dname '/metropolis/' M_.fname '_' type '_draws' num2str(file) ]);
+        temp=load([M_.dname '/metropolis/' M_.fname '_' type '_draws' num2str(file) ]);
     else
-        load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
+        temp=load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
     end
-    NumberOfDraws = rows(pdraws);
-    isdrsaved = columns(pdraws)-1;
+    NumberOfDraws = rows(temp.pdraws);
+    isdrsaved = columns(temp.pdraws)-1;
     for linee = 1:NumberOfDraws
         linea = linea+1;
         if isdrsaved
-            M_=set_parameters_locally(M_,pdraws{linee,1});% Needed to update the covariance matrix of the state innovations.
-            dr = pdraws{linee,2};
+            M_=set_parameters_locally(M_,temp.pdraws{linee,1});% Needed to update the covariance matrix of the state innovations.
+            dr = temp.pdraws{linee,2};
         else
-            M_=set_parameters_locally(M_,pdraws{linee,1});
+            M_=set_parameters_locally(M_,temp.pdraws{linee,1});
             [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
         end
         tmp = th_autocovariances(dr,ivar,M_,options_,nodecomposition);
diff --git a/matlab/dsge_simulated_theoretical_variance_decomposition.m b/matlab/dsge_simulated_theoretical_variance_decomposition.m
index 7e6942c194..cc681eed7c 100644
--- a/matlab/dsge_simulated_theoretical_variance_decomposition.m
+++ b/matlab/dsge_simulated_theoretical_variance_decomposition.m
@@ -85,7 +85,7 @@ NumberOfSavedElementsPerSimulation = nvar*(nexo+1);
 MaXNumberOfDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation/8);
 
 ME_present=0;
-if ~all(M_.H==0)
+if ~all(diag(M_.H)==0)
     if isoctave || matlab_ver_less_than('8.1')
         [observable_pos_requested_vars,index_subset,index_observables]=intersect_stable(ivar,options_.varobs_id);
     else
@@ -130,20 +130,20 @@ linea_ME = 0;
 only_non_stationary_vars=0;
 for file = 1:NumberOfDrawsFiles
     if posterior
-        load([M_.dname '/metropolis/' M_.fname '_' type '_draws' num2str(file) ]);
+        temp=load([M_.dname '/metropolis/' M_.fname '_' type '_draws' num2str(file) ]);
     else
-        load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
+        temp=load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
     end
-    isdrsaved = columns(pdraws)-1;
-    NumberOfDraws = rows(pdraws);
+    isdrsaved = columns(temp.pdraws)-1;
+    NumberOfDraws = rows(temp.pdraws);
     for linee = 1:NumberOfDraws
         linea = linea+1;
         linea_ME = linea_ME+1;
         if isdrsaved
-            M_=set_parameters_locally(M_,pdraws{linee,1});% Needed to update the covariance matrix of the state innovations.
-            dr = pdraws{linee,2};
+            M_=set_parameters_locally(M_,temp.pdraws{linee,1});% Needed to update the covariance matrix of the state innovations.
+            dr = temp.pdraws{linee,2};
         else
-            M_=set_parameters_locally(M_,pdraws{linee,1});
+            M_=set_parameters_locally(M_,temp.pdraws{linee,1});
             [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
         end
         if file==1 && linee==1
@@ -164,6 +164,7 @@ for file = 1:NumberOfDrawsFiles
                 end
             end
             if ME_present
+                M_ = set_measurement_errors(temp.pdraws{linee,1},temp.estim_params_,M_);
                 ME_Variance=diag(M_.H);
                 tmp_ME=NaN(nobs_ME,nexo+1);
                 tmp_ME(:,1:end-1)=tmp{2}(index_subset,:).*repmat(diag(tmp{1}(index_subset,index_subset))./(diag(tmp{1}(index_subset,index_subset))+ME_Variance(index_observables)),1,nexo);
diff --git a/matlab/posterior_analysis.m b/matlab/posterior_analysis.m
index 16088823d2..4d4c50fc3f 100644
--- a/matlab/posterior_analysis.m
+++ b/matlab/posterior_analysis.m
@@ -64,7 +64,7 @@ switch type
     end
     oo_ = variance_decomposition_mc_analysis(SampleSize,'posterior',M_.dname,M_.fname,...
                                              M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
-    if ~all(M_.H==0)
+    if ~all(diag(M_.H)==0)
         if strmatch(arg1,options_.varobs,'exact')
             if isoctave || matlab_ver_less_than('8.1')
                 [observable_name_requested_vars,index_subset,index_observables]=intersect_stable(vartan,options_.varobs);
@@ -72,7 +72,7 @@ switch type
                 [observable_name_requested_vars,index_subset,index_observables]=intersect(vartan,options_.varobs,'stable');
             end
             oo_ = variance_decomposition_ME_mc_analysis(SampleSize,'posterior',M_.dname,M_.fname,...
-                                                        M_.exo_names,arg2,observable_name_requested_vars,arg1,options_.mh_conf_sig,oo_,options_);
+                                                        [M_.exo_names;'ME'],arg2,observable_name_requested_vars,arg1,options_.mh_conf_sig,oo_,options_);
         end
     end
   case 'correlation'
@@ -89,10 +89,10 @@ switch type
     end
     oo_ = conditional_variance_decomposition_mc_analysis(SampleSize,'posterior',M_.dname,M_.fname,...
                                                       arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
-    if ~all(M_.H==0)
+    if ~all(diag(M_.H)==0)
         if strmatch(arg1,options_.varobs,'exact')
             oo_ = conditional_variance_decomposition_ME_mc_analysis(SampleSize,'posterior',M_.dname,M_.fname,...
-                                                              arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
+                                                              arg3,[M_.exo_names;'ME'],arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
         end
     end
   otherwise
diff --git a/matlab/prior_analysis.m b/matlab/prior_analysis.m
index 0c945c3ebd..3a621a9a82 100644
--- a/matlab/prior_analysis.m
+++ b/matlab/prior_analysis.m
@@ -78,7 +78,7 @@ switch type
     end
     oo_ = conditional_variance_decomposition_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
                                                       arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
-    if ~all(M_.H==0)
+    if ~all(diag(M_.H)==0)
         if strmatch(vartan(arg1,:),options_.varobs,'exact')
             oo_ = conditional_variance_decomposition_ME_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
                                                               arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
diff --git a/matlab/prior_sampler.m b/matlab/prior_sampler.m
index 2506f91c30..c09ff8c66d 100644
--- a/matlab/prior_sampler.m
+++ b/matlab/prior_sampler.m
@@ -150,7 +150,7 @@ while iteration < NumberOfSimulations
     end
     if ( file_line_number==TableOfInformations(file_indx_number+1,2) )
         file_indx_number = file_indx_number + 1;
-        save([ PriorDirectoryName '/prior_draws' int2str(file_indx_number) '.mat' ],'pdraws');
+        save([ PriorDirectoryName '/prior_draws' int2str(file_indx_number) '.mat' ],'pdraws','estim_params_');
         if file_indx_number<NumberOfFiles
             if drsave
                 pdraws = cell(TableOfInformations(file_indx_number+1,2),drsave+2);
diff --git a/matlab/selec_posterior_draws.m b/matlab/selec_posterior_draws.m
index 11d9a4dad1..c42e1ced80 100644
--- a/matlab/selec_posterior_draws.m
+++ b/matlab/selec_posterior_draws.m
@@ -122,7 +122,7 @@ if info
             old_mhblck = mhblck;
         end
         clear('x2')
-        save([BaseName '_posterior_draws1.mat'],'pdraws')
+        save([BaseName '_posterior_draws1.mat'],'pdraws','estim_params_')
     else% The posterior draws are saved in xx files.
         NumberOfDrawsPerFile = fix(MAX_mega_bytes/drawsize);
         NumberOfFiles = ceil(SampleSize*drawsize/MAX_mega_bytes);
@@ -149,7 +149,7 @@ if info
             old_mhblck = mhblck;
             if fnum < NumberOfFiles && linee == NumberOfDrawsPerFile
                 linee = 0;
-                save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws')
+                save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws','estim_params_')
                 fnum = fnum+1;
                 if fnum < NumberOfFiles
                     pdraws = cell(NumberOfDrawsPerFile,info);
@@ -158,6 +158,6 @@ if info
                 end
             end
         end
-        save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws')
+        save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws','estim_params_')
     end
 end
\ No newline at end of file
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index 3418770fd5..7081d500b5 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -124,7 +124,7 @@ if ~options_.noprint
         lh = cellofchararraymaxlength(labels)+2;
         dyn_latex_table(M_, options_, my_title, 'covar_ex_shocks', headers, labels, M_.Sigma_e, lh, 10, 6);
     end
-    if ~all(M_.H==0)
+    if ~all(diag(M_.H)==0)
         my_title='MATRIX OF COVARIANCE OF MEASUREMENT ERRORS';
         labels = cellfun(@(x) horzcat('SE_', x), options_.varobs, 'UniformOutput', false);
         headers = vertcat('Variables', labels);
-- 
GitLab