diff --git a/matlab/compute_moments_varendo.m b/matlab/compute_moments_varendo.m
index 21f5201f00984fea5a20d920a0bf6621c95de3cc..d4a35758b6330b15936528e88785d2322887c9b4 100644
--- a/matlab/compute_moments_varendo.m
+++ b/matlab/compute_moments_varendo.m
@@ -42,6 +42,9 @@ if strcmpi(type,'posterior')
     if nargin==4
         var_list_ = options_.varobs;
     end
+    if isfield(oo_,'PosteriorTheoreticalMoments')
+        oo_=rmfield(oo_,'PosteriorTheoreticalMoments');
+    end
 elseif strcmpi(type,'prior')
     posterior = 0;
     if nargin==4
@@ -50,6 +53,9 @@ elseif strcmpi(type,'prior')
             options_.prior_analysis_var_list = options_.varobs;
         end
     end
+    if isfield(oo_,'PriorTheoreticalMoments')
+        oo_=rmfield(oo_,'PriorTheoreticalMoments');
+    end
 else
     error('compute_moments_varendo:: Unknown type!')
 end
@@ -144,7 +150,7 @@ if M_.exo_nbr > 1
         skipline();
     end
     skipline();
-    if ~all(M_.H==0)
+    if ~all(diag(M_.H)==0)
         if isoctave
             [observable_name_requested_vars, varlist_pos] = intersect_stable(var_list_, options_.varobs);
         else
@@ -231,7 +237,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 45fe803bdc89f149a35931aa521b4fc20636bf94..f249a7c7a9e17b7de25b061c2ae30bdbe3cbaa35 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
         [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 6fd28a0d654ffe3a7ac19a717b821d96c86b9fa9..1f71d9c9549b27d19c66f894e7a58e1ae85e75c1 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
         [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 4195b58a91871168590a346df97435be8ac71a4d..22ced1fe5c37d79defee7d76cb1133054750030f 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
         [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 fb00d5ddcd537074f918a8cc63c565a5b79da4d4..c32aa9db4226eb65bf8dd821180867ceb76017c7 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 0be948bc3c2f7f01601853278939e3af91058637..bbd4702751335d3d918386a26aa4b7c3a2f17908 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 71febcf473cb55e8399541910b453c2ff94852fd..b93aa246ed7f5370f8dc7bec05aef5eb4be4a41a 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
         [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 65735025c20b90addf0c4026038463dfaf03efa5..03bff87b8733cae15839ba2351a97d1c797db88a 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
                 [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 0c945c3ebd6d603fd8927bfe598a979bb8731116..3a621a9a8237dd88486f392926c8e0e7128de91f 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 2506f91c306aef73867c2b982ebd8736b0e3615a..c09ff8c66de7a427480afed64aa86efcd24ef48a 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 11d9a4dad16e6388f81134ba135670f8e91839ce..c42e1ced80fdc769dfe10a731cc0e2a45483c7e9 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/set_measurement_errors.m b/matlab/set_measurement_errors.m
new file mode 100644
index 0000000000000000000000000000000000000000..c51d7c4772db0fa1ce37e4610bb183d83ccea85c
--- /dev/null
+++ b/matlab/set_measurement_errors.m
@@ -0,0 +1,71 @@
+function M_ = set_measurement_errors(xparam1,estim_params_,M_)
+% function M_=set_measurement_errors(xparam1,estim_params_,M_)
+% Sets parameters value (except measurement errors)
+% This is called for computations such as IRF and forecast
+% when measurement errors aren't taken into account; in contrast to
+% set_parameters.m, the global M_-structure is not altered
+%
+% INPUTS
+%    xparam1:   vector of parameters to be estimated (initial values)
+%    M_:        Dynare model-structure
+%
+% OUTPUTS
+%    M_:        Dynare model-structure
+%
+% SPECIAL REQUIREMENTS
+%    none
+
+% Copyright (C) 2017 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+H = M_.H;
+Correlation_matrix_ME = M_.Correlation_matrix_ME;
+
+% setting measument error variance; on the diagonal of Covariance matrix; used later
+% for updating covariances
+offset = estim_params_.nvx;
+
+if estim_params_.nvn
+    for i=1:estim_params_.nvn
+        k = estim_params_.nvn_observable_correspondence(i,1);
+        H(k,k) = xparam1(i+offset)^2;
+    end
+end
+
+% update offset
+offset = estim_params_.nvx+estim_params_.nvn+estim_params_.ncx;
+
+% setting measurement error covariances
+if estim_params_.ncn
+    corrn_observable_correspondence = estim_params_.corrn_observable_correspondence;
+    for i=1:estim_params_.ncn
+        k1 = corrn_observable_correspondence(i,1);
+        k2 = corrn_observable_correspondence(i,2);
+        Correlation_matrix_ME(k1,k2) = xparam1(i+offset);
+        Correlation_matrix_ME(k2,k1) = Correlation_matrix_ME(k1,k2);
+    end
+end
+%build covariance matrix from correlation matrix and variances already on
+%diagonal
+H = diag(sqrt(diag(H)))*Correlation_matrix_ME*diag(sqrt(diag(H)));
+%if calibrated covariances, set them now to their stored value
+if isfield(estim_params_,'calibrated_covariances_ME')
+    H(estim_params_.calibrated_covariances_ME.position)=estim_params_.calibrated_covariances_ME.cov_value;
+end
+
+M_.H = H;
+M_.Correlation_matrix_ME=Correlation_matrix_ME;
diff --git a/matlab/set_parameters_locally.m b/matlab/set_parameters_locally.m
index 492024daa8e503ecfd913624378656662e275013..008abc3a0498192daed277ba0d586a5daf6612b3 100644
--- a/matlab/set_parameters_locally.m
+++ b/matlab/set_parameters_locally.m
@@ -1,6 +1,6 @@
 function M_=set_parameters_locally(M_,xparam1)
 
-% function M_out=set_parameters(M_,xparam1)
+% function M_=set_parameters_locally(M_,xparam1)
 % Sets parameters value (except measurement errors)
 % This is called for computations such as IRF and forecast
 % when measurement errors aren't taken into account; in contrast to
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index ac0908de6b8e102249d63dcd30476c8d7f4ea0e2..2170f0027586e45bec6599096f403bcef53f2094 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);
diff --git a/tests/estimation/fs2000_calibrated_covariance.mod b/tests/estimation/fs2000_calibrated_covariance.mod
index 79d309786857d00238fc1fb57a88e1cdf4972b20..a2acc3387567c23ed748063b5981b9afe5bb7b2c 100644
--- a/tests/estimation/fs2000_calibrated_covariance.mod
+++ b/tests/estimation/fs2000_calibrated_covariance.mod
@@ -80,7 +80,7 @@ corr e_a, e_m, 0.5;
 stderr gp_obs, 0.5;
 end;
 
-estimation(order=1,datafile=fsdat_simul,nobs=192, loglinear, mh_replic=2002, mh_nblocks=1, mh_jscale=0.8);
+estimation(order=1,datafile=fsdat_simul,nobs=192, loglinear, mh_replic=0, mh_nblocks=1, mh_jscale=0.8,moments_varendo,consider_all_endogenous);
 
 if isequal(M_.Sigma_e(2,1),5e-5) || isequal(M_.Sigma_e(1,2),5e-5)
     error('Problem in overriding calibrated covariance of structural shocks by estimated correlation')
@@ -88,3 +88,5 @@ end
 if ~isequal(M_.H(2,1),5e-5) || ~isequal(M_.H(1,2),5e-5)
     error('Problem in setting calibrated covariance of measurement errors')
 end
+
+stoch_simul(order=1,periods=1000);
diff --git a/tests/kalman_filter_smoother/check_variable_dimensions/fs2000_ML.mod b/tests/kalman_filter_smoother/check_variable_dimensions/fs2000_ML.mod
index b333d7ff0773dd523b5e33ef114999a94e353251..19d5c5c1a24da0f247a95628c7450bc1fa4fa02e 100644
--- a/tests/kalman_filter_smoother/check_variable_dimensions/fs2000_ML.mod
+++ b/tests/kalman_filter_smoother/check_variable_dimensions/fs2000_ML.mod
@@ -116,7 +116,7 @@ corr e_m, e_a, 0;
 stderr gp_obs, 0.01;
 end;
 options_.prior_trunc=0;
-estimation(order=1,datafile='../fsdat_simul', nobs=192, loglinear, forecast=8,smoother,filter_covariance,filtered_vars,filter_step_ahead=[1,2,4],filter_decomposition,selected_variables_only)  m P c e W R k d y gy_obs;
+estimation(order=1,datafile='../fsdat_simul', nobs=192, loglinear, moments_varendo,conditional_variance_decomposition=[1,3],forecast=8,smoother,filter_covariance,filtered_vars,filter_step_ahead=[1,2,4],filter_decomposition,selected_variables_only)  m P c e W R k d y gy_obs gp_obs;
 
 
 if size(oo_.FilteredVariablesKStepAhead,3)~=(options_.nobs+max(options_.filter_step_ahead)) || ...
diff --git a/tests/moments/example1_one_sided_hp_test.mod b/tests/moments/example1_one_sided_hp_test.mod
index 1e61f23d26ef407864f1dbfdecfc3bc95f944120..c87763ff06c2a01f9a1572c2945cd56f781d0345 100644
--- a/tests/moments/example1_one_sided_hp_test.mod
+++ b/tests/moments/example1_one_sided_hp_test.mod
@@ -68,3 +68,14 @@ end;
 steady(solve_algo=4,maxit=1000);
 
 stoch_simul(order=1,nofunctions,one_sided_hp_filter=1600,irf=0,periods=5000,filtered_theoretical_moments_grid=8192);
+
+varobs k c y;
+
+shocks;
+var e; stderr 0.009;
+var u; stderr 0.009;
+var e, u = phi*0.009*0.009;
+var c; stderr 0.01;
+end;
+
+stoch_simul(order=1,nofunctions,one_sided_hp_filter=1600,irf=0,periods=5000,filtered_theoretical_moments_grid=8192);
diff --git a/tests/moments/fs2000_post_moments.mod b/tests/moments/fs2000_post_moments.mod
index 3ea15dcef0946954f34f1a85f03927a40fa13c09..fce555f413c5c50a86e4267c5a0ca679ca4509f6 100644
--- a/tests/moments/fs2000_post_moments.mod
+++ b/tests/moments/fs2000_post_moments.mod
@@ -169,7 +169,7 @@ end
 
 conditional_variance_decomposition=mean(conditional_variance_decomposition,4);
 nvars=M_.orig_endo_nbr;
-horizon_size=size(conditional_variance_decomposition,3);
+horizon_size=size(conditional_variance_decomposition,2);
 for var_iter_1=1:nvars
     for shock_iter=1:M_.exo_nbr
         for horizon_iter=1:horizon_size
@@ -180,6 +180,105 @@ for var_iter_1=1:nvars
     end
 end
 
+// case with measurement error
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.100;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+stderr gp_obs, inv_gamma_pdf, 0.003, inf;
+end;
+
+estimation(order=1,mode_compute=5, datafile='../fs2000/fsdat_simul.m', nobs=192, loglinear, mh_replic=20, mh_nblocks=1, mh_jscale=0.8,moments_varendo,
+conditional_variance_decomposition=[2,2000],consider_all_endogenous,sub_draws=2);
+
+stoch_simul(order=1,conditional_variance_decomposition=[2,2000],noprint,nograph);
+par=load([M_.fname filesep 'metropolis' filesep M_.fname '_posterior_draws1']);
+
+for par_iter=1:size(par.pdraws,1)
+   M_=set_parameters_locally(M_,par.pdraws{par_iter,1});
+   [info, oo_, options_, M_]=stoch_simul(M_, options_, oo_, var_list_);
+   correlation(:,:,par_iter)=cell2mat(oo_.autocorr);
+   covariance(:,:,par_iter)=oo_.var;
+   conditional_variance_decomposition(:,:,:,par_iter)=oo_.conditional_variance_decomposition;
+   conditional_variance_decomposition_ME(:,:,:,par_iter)=oo_.conditional_variance_decomposition_ME;
+   variance_decomposition(:,:,par_iter)=oo_.variance_decomposition;
+   variance_decomposition_ME(:,:,par_iter)=oo_.variance_decomposition_ME;
+   [~,obs_order]=sort(options_.varobs_id);
+end
+
+correlation=mean(correlation,3);
+nvars=M_.orig_endo_nbr;
+for var_iter_1=1:nvars
+    for var_iter_2=1:nvars
+        if max(abs(correlation(var_iter_1,var_iter_2:nvars:end)'-oo_.PosteriorTheoreticalMoments.dsge.correlation.Mean.(M_.endo_names{var_iter_1}).(M_.endo_names{var_iter_2})))>1e-8
+            error('Correlations do not match')
+        end
+    end
+end
+
+covariance=mean(covariance,3);
+nvars=M_.orig_endo_nbr;
+for var_iter_1=1:nvars
+    for var_iter_2=var_iter_1:nvars
+        if max(abs(covariance(var_iter_1,var_iter_2)-oo_.PosteriorTheoreticalMoments.dsge.covariance.Mean.(M_.endo_names{var_iter_1}).(M_.endo_names{var_iter_2})))>1e-8
+            error('Covariances do not match')
+        end
+    end
+end
+
+variance_decomposition=mean(variance_decomposition,3);
+nvars=M_.orig_endo_nbr;
+for var_iter_1=1:nvars
+    for shock_iter=1:M_.exo_nbr
+        if max(abs(variance_decomposition(var_iter_1,shock_iter)/100-oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.Mean.(M_.endo_names{var_iter_1}).(M_.exo_names{shock_iter})))>1e-8
+            error('Variance decomposition does not match')
+        end
+    end
+end
+
+variance_decomposition_ME=mean(variance_decomposition_ME,3);
+nvars=length(options_.varobs);
+for var_iter_1=1:nvars
+    for shock_iter=1:M_.exo_nbr
+        if max(abs(variance_decomposition_ME(obs_order(var_iter_1),shock_iter)/100-oo_.PosteriorTheoreticalMoments.dsge.VarianceDecompositionME.Mean.(options_.varobs{var_iter_1}).(M_.exo_names{shock_iter})))>1e-8
+            error('Variance decomposition does not match')
+        end
+    end
+end
+
+conditional_variance_decomposition=mean(conditional_variance_decomposition,4);
+nvars=M_.orig_endo_nbr;
+horizon_size=size(conditional_variance_decomposition,2);
+for var_iter_1=1:nvars
+    for shock_iter=1:M_.exo_nbr
+        for horizon_iter=1:horizon_size
+            if max(abs(conditional_variance_decomposition(var_iter_1,horizon_iter,shock_iter)-oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecomposition.Mean.(M_.endo_names{var_iter_1}).(M_.exo_names{shock_iter})(horizon_iter)))>1e-8
+                error('Conditional Variance decomposition does not match')
+            end
+        end
+    end
+end
+
+conditional_variance_decomposition_ME=mean(conditional_variance_decomposition_ME,4);
+exo_names=[M_.exo_names;'ME'];
+nvars=length(options_.varobs);
+horizon_size=size(conditional_variance_decomposition_ME,2);
+for var_iter_1=1:nvars
+    for shock_iter=1:M_.exo_nbr+1
+        for horizon_iter=1:horizon_size
+            if max(abs(conditional_variance_decomposition_ME(obs_order(var_iter_1),horizon_iter,shock_iter)-oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecompositionME.Mean.(options_.varobs{var_iter_1}).(exo_names{shock_iter})(horizon_iter)))>1e-8
+                error('Conditional Variance decomposition does not match')
+            end
+        end
+    end
+end
+
 /*
  * The following lines were used to generate the data file. If you want to
  * generate another random data file, comment the "estimation" line and uncomment