From 6e06acc7f4fd8fed9f8dfc647a021770e952ed15 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer 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 --- 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 21f5201f0..052d62217 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 [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 45fe803bd..f249a7c7a 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 6fd28a0d6..1f71d9c95 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 4195b58a9..22ced1fe5 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 fb00d5ddc..c32aa9db4 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 0be948bc3..bbd470275 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 71febcf47..b93aa246e 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 65735025c..03bff87b8 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 0c945c3eb..3a621a9a8 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 2506f91c3..c09ff8c66 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