From 1ef1f9c75d4adece4fac9eecc192356e62b6f582 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 24 Aug 2015 22:27:10 +0200 Subject: [PATCH 1/3] Transform hard-coded kernel density options to real options --- matlab/global_initialization.m | 4 ++++ matlab/posterior_moments.m | 15 +++++++++++---- 2 files changed, 15 insertions(+), 4 deletions(-) diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 3f65232c9..14550cc85 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -487,6 +487,10 @@ options_.selected_variables_only = 0; options_.contemporaneous_correlation = 0; options_.initialize_estimated_parameters_with_the_prior_mode = 0; options_.estimation_dll = 0; +options_.estimation.posterior_kernel_density.indicator = 0; +options_.estimation.posterior_kernel_density.gridpoints = 2^9; +options_.estimation.posterior_kernel_density.bandwidth = 0; % Rule of thumb optimal bandwidth parameter. +options_.estimation.posterior_kernel_density.kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton. % Misc options_.conf_sig = 0.6; oo_.exo_simul = []; diff --git a/matlab/posterior_moments.m b/matlab/posterior_moments.m index 7e122cf46..beb1b0a41 100644 --- a/matlab/posterior_moments.m +++ b/matlab/posterior_moments.m @@ -1,4 +1,4 @@ -function [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(xx,info,mh_conf_sig) +function [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(xx,info,mh_conf_sig,kernel_options) % Computes posterior mean, median, variance, HPD interval, deciles, and density from posterior draws. % % INPUTS @@ -38,6 +38,16 @@ function [post_mean, post_median, post_var, hpd_interval, post_deciles, density] % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . +if nargin<4 + number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two. + bandwidth = 0; % Rule of thumb optimal bandwidth parameter. + kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton. +else + number_of_grid_points = kernel_options.gridpoints; + bandwidth = kernel_options.bandwidth; + kernel_function = kernel_options.kernel_function; +end + xx = xx(:); xx = sort(xx); @@ -76,10 +86,7 @@ else end density = []; if info - number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two. if post_var > 1e-12 - bandwidth = 0; % Rule of thumb optimal bandwidth parameter. - kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton. optimal_bandwidth = mh_optimal_bandwidth(xx,number_of_draws,bandwidth,kernel_function); [density(:,1),density(:,2)] = kernel_density_estimate(xx,number_of_grid_points,... number_of_draws,optimal_bandwidth,kernel_function); From b7cbb563d67491869f28c8c0069ac13c2bda45f3 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Tue, 25 Aug 2015 09:22:40 +0200 Subject: [PATCH 2/3] Allow suppressing density of smoother and forecast objects --- doc/dynare.texi | 13 +++- ...ional_variance_decomposition_mc_analysis.m | 21 ++++-- matlab/correlation_mc_analysis.m | 37 ++++++---- matlab/covariance_mc_analysis.m | 26 +++++-- matlab/global_initialization.m | 8 +-- matlab/pm3.m | 68 ++++++++++++------- matlab/posterior_analysis.m | 4 +- matlab/prior_analysis.m | 4 +- matlab/variance_decomposition_mc_analysis.m | 22 ++++-- preprocessor/DynareBison.yy | 8 +-- preprocessor/DynareFlex.ll | 2 +- tests/estimation/fs2000.mod | 2 +- 12 files changed, 142 insertions(+), 73 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index 4f2d3f4ee..e60b530fb 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -5436,6 +5436,9 @@ controlled by the @code{ar} option. @xref{contemporaneous_correlation}. Results are stored in @code{oo_.PosteriorTheoreticalMoments}. Note that the @code{nocorr}-option has no effect. +@item no_posterior_kernel_density +Shuts off the computation of the kernel density estimator for the posterior objects (@pxref{density-field}). + @item conditional_variance_decomposition = @var{INTEGER} See below. @@ -5814,8 +5817,10 @@ Variance of the posterior distribution Deciles of the distribution. @item density -Non parametric estimate of the posterior density. First and second -columns are respectively abscissa and ordinate coordinates. +@anchor{density-field} +Non parametric estimate of the posterior density following the approach outlined in @cite{Skoeld and Roberts (2003)}. First and second +columns are respectively abscissa and ordinate coordinates. + @end table @item ESTIMATED_OBJECT @@ -13124,6 +13129,10 @@ Sims, Christopher A., Daniel F. Waggoner and Tao Zha (2008): ``Methods for inference in large multiple-equation Markov-switching models,'' @i{Journal of Econometrics}, 146, 255--274 +@item +Skoeld, Martin and Gareth O. Roberts (2003): ``Density Estimation for the +Metropolis-Hastings Algorithm,'' @i{Scandinavian Journal of Statistics}, 30, 699--718 + @item Smets, Frank and Rafael Wouters (2003): ``An Estimated Dynamic Stochastic General Equilibrium Model of the Euro Area,'' @i{Journal of diff --git a/matlab/conditional_variance_decomposition_mc_analysis.m b/matlab/conditional_variance_decomposition_mc_analysis.m index 22273f2be..2c724083a 100644 --- a/matlab/conditional_variance_decomposition_mc_analysis.m +++ b/matlab/conditional_variance_decomposition_mc_analysis.m @@ -1,5 +1,5 @@ function oo_ = ... - conditional_variance_decomposition_mc_analysis(NumberOfSimulations, type, dname, fname, Steps, exonames, exo, var_list, endogenous_variable_index, mh_conf_sig, oo_) + conditional_variance_decomposition_mc_analysis(NumberOfSimulations, type, dname, fname, Steps, exonames, exo, var_list, endogenous_variable_index, mh_conf_sig, oo_,options_) % This function analyses the (posterior or prior) distribution of the % endogenous variables' conditional variance decomposition. % @@ -91,19 +91,26 @@ p_mean = NaN(1,length(Steps)); p_median = NaN(1,length(Steps)); p_variance = NaN(1,length(Steps)); p_deciles = NaN(9,length(Steps)); -p_density = NaN(2^9,2,length(Steps)); +if options_.estimation.moments_posterior_density.indicator + p_density = NaN(2^9,2,length(Steps)); +end p_hpdinf = NaN(1,length(Steps)); p_hpdsup = NaN(1,length(Steps)); for i=1:length(Steps) - [pp_mean, pp_median, pp_var, hpd_interval, pp_deciles, pp_density] = ... - posterior_moments(tmp(:,i),1,mh_conf_sig); + if options_.estimation.moments_posterior_density.indicator + [pp_mean, pp_median, pp_var, hpd_interval, pp_deciles, pp_density] = ... + posterior_moments(tmp(:,i),1,mh_conf_sig); + p_density(:,:,i) = pp_density; + else + [pp_mean, pp_median, pp_var, hpd_interval, pp_deciles] = ... + posterior_moments(tmp(:,i),0,mh_conf_sig); + end p_mean(i) = pp_mean; p_median(i) = pp_median; p_variance(i) = pp_var; p_deciles(:,i) = pp_deciles; p_hpdinf(i) = hpd_interval(1); p_hpdsup(i) = hpd_interval(2); - p_density(:,:,i) = pp_density; end eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.Steps = Steps;']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.Mean.' name ' = p_mean;']); @@ -112,4 +119,6 @@ eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.Vari eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.HPDinf.' name ' = p_hpdinf;']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.HPDsup.' name ' = p_hpdsup;']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.deciles.' name ' = p_deciles;']); -eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.density.' name ' = p_density;']); \ No newline at end of file +if options_.estimation.moments_posterior_density.indicator + eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.density.' name ' = p_density;']); +end \ No newline at end of file diff --git a/matlab/correlation_mc_analysis.m b/matlab/correlation_mc_analysis.m index 4654b50bd..28a5c1af7 100644 --- a/matlab/correlation_mc_analysis.m +++ b/matlab/correlation_mc_analysis.m @@ -67,19 +67,19 @@ if isfield(oo_,[TYPE 'TheoreticalMoments']) end end else - oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_); + oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_,options_); end else - oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_); + oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_,options_); end else - oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_); + oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_,options_); end else - oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_); + oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_,options_); end else - oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_); + oo_ = initialize_output_structure(var1,var2,nar,TYPE,oo_,options_); end ListOfFiles = dir([ PATH fname '_' TYPE 'Correlations*.mat']); i1 = 1; tmp = zeros(SampleSize,1); @@ -91,8 +91,13 @@ for file = 1:length(ListOfFiles) end name = [ var1 '.' var2 ]; if ~isconst(tmp) - [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ... - posterior_moments(tmp,1,mh_conf_sig); + if options_.estimation.moments_posterior_density.indicator + [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ... + posterior_moments(tmp,1,mh_conf_sig); + else + [p_mean, p_median, p_var, hpd_interval, p_deciles] = ... + posterior_moments(tmp,0,mh_conf_sig); + end if isfield(oo_,[ TYPE 'TheoreticalMoments']) eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments;']) if isfield(temporary_structure,'dsge') @@ -104,7 +109,9 @@ if ~isconst(tmp) oo_ = fill_output_structure(var1,var2,TYPE,oo_,'HPDinf',nar,hpd_interval(1)); oo_ = fill_output_structure(var1,var2,TYPE,oo_,'HPDsup',nar,hpd_interval(2)); oo_ = fill_output_structure(var1,var2,TYPE,oo_,'deciles',nar,p_deciles); - oo_ = fill_output_structure(var1,var2,TYPE,oo_,'density',nar,density); + if options_.estimation.moments_posterior_density.indicator + oo_ = fill_output_structure(var1,var2,TYPE,oo_,'density',nar,density); + end end end end @@ -120,13 +127,15 @@ else oo_ = fill_output_structure(var1,var2,TYPE,oo_,'HPDinf',nar,NaN); oo_ = fill_output_structure(var1,var2,TYPE,oo_,'HPDsup',nar,NaN); oo_ = fill_output_structure(var1,var2,TYPE,oo_,'deciles',nar,NaN); - oo_ = fill_output_structure(var1,var2,TYPE,oo_,'density',nar,NaN); + if options_.estimation.moments_posterior_density.indicator + oo_ = fill_output_structure(var1,var2,TYPE,oo_,'density',nar,NaN); + end end end end end -function oo_ = initialize_output_structure(var1,var2,nar,type,oo_) +function oo_ = initialize_output_structure(var1,var2,nar,type,oo_,options_) name = [ var1 '.' var2 ]; eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.Mean.' name ' = NaN(' int2str(nar) ',1);']); eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.Median.' name ' = NaN(' int2str(nar) ',1);']); @@ -134,9 +143,13 @@ eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.Variance.' name ' = NaN(' eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.HPDinf.' name ' = NaN(' int2str(nar) ',1);']); eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.HPDsup.' name ' = NaN(' int2str(nar) ',1);']); eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.deciles.' name ' = cell(' int2str(nar) ',1);']); -eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.density.' name ' = cell(' int2str(nar) ',1);']); +if options_.estimation.moments_posterior_density.indicator + eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.density.' name ' = cell(' int2str(nar) ',1);']); +end for i=1:nar - eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.density.' name '(' int2str(i) ',1) = {NaN};']); + if options_.estimation.moments_posterior_density.indicator + eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.density.' name '(' int2str(i) ',1) = {NaN};']); + end eval(['oo_.' type 'TheoreticalMoments.dsge.correlation.deciles.' name '(' int2str(i) ',1) = {NaN};']); end diff --git a/matlab/covariance_mc_analysis.m b/matlab/covariance_mc_analysis.m index 59cc635f2..991c9f1ac 100644 --- a/matlab/covariance_mc_analysis.m +++ b/matlab/covariance_mc_analysis.m @@ -106,15 +106,20 @@ for file = 1:length(ListOfFiles) end name = [var1 '.' var2]; if ~isconst(tmp) - [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ... - posterior_moments(tmp,1,mh_conf_sig); + if options_.estimation.moments_posterior_density.indicator + [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ... + posterior_moments(tmp,1,mh_conf_sig); + eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.density.' name ' = density;']); + else + [p_mean, p_median, p_var, hpd_interval, p_deciles] = ... + posterior_moments(tmp,0,mh_conf_sig); + end eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.Mean.' name ' = p_mean;']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.Median.' name ' = p_median;']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.Variance.' name ' = p_var;']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.HPDinf.' name ' = hpd_interval(1);']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.HPDsup.' name ' = hpd_interval(2);']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.deciles.' name ' = p_deciles;']); - eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.density.' name ' = density;']); else eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.Mean.' name ' = NaN;']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.Median.' name ' = NaN;']); @@ -122,17 +127,24 @@ else eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.HPDinf.' name ' = NaN;']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.HPDsup.' name ' = NaN;']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.deciles.' name ' = NaN;']); - eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.density.' name ' = NaN;']); + if options_.estimation.moments_posterior_density.indicator + eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.density.' name ' = NaN;']); + end end if options_.contemporaneous_correlation - [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ... - posterior_moments(tmp_corr_mat,1,mh_conf_sig); + if options_.estimation.moments_posterior_density.indicator + [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ... + posterior_moments(tmp_corr_mat,1,mh_conf_sig); + eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.density.' name ' = density;']); + else + [p_mean, p_median, p_var, hpd_interval, p_deciles] = ... + posterior_moments(tmp_corr_mat,0,mh_conf_sig); + end eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.Mean.' name ' = p_mean;']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.Median.' name ' = p_median;']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.Variance.' name ' = p_var;']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.HPDinf.' name ' = hpd_interval(1);']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.HPDsup.' name ' = hpd_interval(2);']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.deciles.' name ' = p_deciles;']); - eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.density.' name ' = density;']); end diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 14550cc85..41142f29b 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -487,10 +487,10 @@ options_.selected_variables_only = 0; options_.contemporaneous_correlation = 0; options_.initialize_estimated_parameters_with_the_prior_mode = 0; options_.estimation_dll = 0; -options_.estimation.posterior_kernel_density.indicator = 0; -options_.estimation.posterior_kernel_density.gridpoints = 2^9; -options_.estimation.posterior_kernel_density.bandwidth = 0; % Rule of thumb optimal bandwidth parameter. -options_.estimation.posterior_kernel_density.kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton. +options_.estimation.moments_posterior_density.indicator = 1; +options_.estimation.moments_posterior_density.gridpoints = 2^9; +options_.estimation.moments_posterior_density.bandwidth = 0; % Rule of thumb optimal bandwidth parameter. +options_.estimation.moments_posterior_density.kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton. % Misc options_.conf_sig = 0.6; oo_.exo_simul = []; diff --git a/matlab/pm3.m b/matlab/pm3.m index 21e46792f..8e98004cb 100644 --- a/matlab/pm3.m +++ b/matlab/pm3.m @@ -79,6 +79,9 @@ Median = zeros(n2,nvar); Var = zeros(n2,nvar); Distrib = zeros(9,n2,nvar); HPD = zeros(2,n2,nvar); +if options_.estimation.moments_posterior_density.indicator + Density = zeros(options_.estimation.moments_posterior_density.gridpoints,2,n2,nvar); +end fprintf(['Estimation::mcmc: ' tit1 '\n']); stock1 = zeros(n1,n2,B); k = 0; @@ -112,18 +115,33 @@ if filter_step_ahead_indicator Var_filter_step_ahead = zeros(filter_steps,nvar,n2); Distrib_filter_step_ahead = zeros(9,filter_steps,nvar,n2); HPD_filter_step_ahead = zeros(2,filter_steps,nvar,n2); + if options_.estimation.moments_posterior_density.indicator + Density_filter_step_ahead = zeros(options_.estimation.moments_posterior_density.gridpoints,2,filter_steps,nvar,n2); + end end tmp =zeros(B,1); for i = 1:nvar for j = 1:n2 - [Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i)] = ... - posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),0,options_.mh_conf_sig); + if options_.estimation.moments_posterior_density.indicator + [Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i),Density(:,:,j,i)] = ... + posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density); + else + [Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i)] = ... + posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),0,options_.mh_conf_sig); + end if filter_step_ahead_indicator - for K_step = 1:length(options_.filter_step_ahead) - [Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j)] = ... - posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),0,options_.mh_conf_sig); - end + if options_.estimation.moments_posterior_density.indicator + for K_step = 1:length(options_.filter_step_ahead) + [Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j),Density_filter_step_ahead(:,:,K_step,i,j) ] = ... + posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density); + end + else + for K_step = 1:length(options_.filter_step_ahead) + [Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j)] = ... + posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),0,options_.mh_conf_sig); + end + end end end end @@ -142,21 +160,27 @@ end for i = 1:nvar name = deblank(names1(SelecVariables(i),:)); - eval(['oo_.' name3 '.Mean.' name ' = Mean(:,i);']); - eval(['oo_.' name3 '.Median.' name ' = Median(:,i);']); - eval(['oo_.' name3 '.Var.' name ' = Var(:,i);']); - eval(['oo_.' name3 '.deciles.' name ' = Distrib(:,:,i);']); - eval(['oo_.' name3 '.HPDinf.' name ' = HPD(1,:,i)'';']); - eval(['oo_.' name3 '.HPDsup.' name ' = HPD(2,:,i)'';']); + oo_.(name3).Mean.(name) = Mean(:,i); + oo_.(name3).Median.(name) = Median(:,i); + oo_.(name3).Var.(name) = Var(:,i); + oo_.(name3).deciles.(name) = Distrib(:,:,i); + oo_.(name3).HPDinf.(name) = HPD(1,:,i)'; + oo_.(name3).HPDsup.(name) = HPD(2,:,i)'; + if options_.estimation.moments_posterior_density.indicator + oo_.(name3).density.(name) = Density(:,:,:,i); + end if filter_step_ahead_indicator for K_step = 1:length(options_.filter_step_ahead) name4=['Filtered_Variables_',num2str(K_step),'_step_ahead']; - eval(['oo_.' name4 '.Mean.' name ' = squeeze(Mean_filter_step_ahead(K_step,i,:));']); - eval(['oo_.' name4 '.Median.' name ' = squeeze(Median_filter_step_ahead(K_step,i,:));']); - eval(['oo_.' name4 '.Var.' name ' = squeeze(Var_filter_step_ahead(K_step,i,:));']); - eval(['oo_.' name4 '.deciles.' name ' = squeeze(Distrib_filter_step_ahead(:,K_step,i,:));']); - eval(['oo_.' name4 '.HPDinf.' name ' = squeeze(HPD_filter_step_ahead(1,K_step,i,:));']); - eval(['oo_.' name4 '.HPDsup.' name ' = squeeze(HPD_filter_step_ahead(2,K_step,i,:));']); + oo_.(name4).Mean.(name) = squeeze(Mean_filter_step_ahead(K_step,i,:)); + oo_.(name4).Median.(name) = squeeze(Median_filter_step_ahead(K_step,i,:)); + oo_.(name4).Var.(name) = squeeze(Var_filter_step_ahead(K_step,i,:)); + oo_.(name4).deciles.(name) = squeeze(Distrib_filter_step_ahead(:,K_step,i,:)); + oo_.(name4).HPDinf.(name) = squeeze(HPD_filter_step_ahead(1,K_step,i,:)); + oo_.(name4).HPDsup.(name) = squeeze(HPD_filter_step_ahead(2,K_step,i,:)); + if options_.estimation.moments_posterior_density.indicator + oo_.(name4).density.(name) = squeeze(Density_filter_step_ahead(:,:,K_step,i,:)); + end end end end @@ -267,10 +291,4 @@ if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) fclose(fidTeX); end -fprintf(['Estimation::mcmc: ' tit1 ', done!\n']); - - - - - - +fprintf(['Estimation::mcmc: ' tit1 ', done!\n']); \ No newline at end of file diff --git a/matlab/posterior_analysis.m b/matlab/posterior_analysis.m index 5209cb3f8..faa117c03 100644 --- a/matlab/posterior_analysis.m +++ b/matlab/posterior_analysis.m @@ -62,7 +62,7 @@ switch type dsge_simulated_theoretical_variance_decomposition(SampleSize,M_,options_,oo_,'posterior'); end oo_ = variance_decomposition_mc_analysis(SampleSize,'posterior',M_.dname,M_.fname,... - M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_); + M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_); case 'correlation' if nargin==narg1 [nvar,vartan,NumberOfFiles] = ... @@ -76,7 +76,7 @@ switch type dsge_simulated_theoretical_conditional_variance_decomposition(SampleSize,arg3,M_,options_,oo_,'posterior'); 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_); + arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_); otherwise disp('Not yet implemented') end \ No newline at end of file diff --git a/matlab/prior_analysis.m b/matlab/prior_analysis.m index b942a509f..fd79414e5 100644 --- a/matlab/prior_analysis.m +++ b/matlab/prior_analysis.m @@ -63,7 +63,7 @@ switch type dsge_simulated_theoretical_variance_decomposition(SampleSize,M_,options_,oo_,'prior'); end oo_ = variance_decomposition_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,... - M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_); + M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_); case 'correlation' if nargin==narg1 [nvar,vartan,NumberOfFiles] = ... @@ -77,7 +77,7 @@ switch type dsge_simulated_theoretical_conditional_variance_decomposition(SampleSize,arg3,M_,options_,oo_,'prior'); 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_); + arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_); otherwise disp('Not yet implemented') end \ No newline at end of file diff --git a/matlab/variance_decomposition_mc_analysis.m b/matlab/variance_decomposition_mc_analysis.m index 798b28afe..7d426bd40 100644 --- a/matlab/variance_decomposition_mc_analysis.m +++ b/matlab/variance_decomposition_mc_analysis.m @@ -1,14 +1,14 @@ -function oo_ = variance_decomposition_mc_analysis(NumberOfSimulations,type,dname,fname,exonames,exo,vartan,var,mh_conf_sig,oo_) +function oo_ = variance_decomposition_mc_analysis(NumberOfSimulations,type,dname,fname,exonames,exo,vartan,var,mh_conf_sig,oo_,options_) % function oo_ = variance_decomposition_mc_analysis(NumberOfSimulations,type,dname,fname,exonames,exo,vartan,var,mh_conf_sig,oo_) % This function analyses the (posterior or prior) distribution of the % endogenous variables' variance decomposition. -% +% % INPUTS % NumberOfSimulations [integer] scalar, number of simulations. % type [string] 'prior' or 'posterior' % dname [string] directory name where to save % fname [string] name of the mod-file -% exonames [string] (n_exo*char_length) character array with names of exogenous variables +% exonames [string] (n_exo*char_length) character array with names of exogenous variables % exo [string] name of current exogenous % variable % vartan [string] (n_endo*char_length) character array with name @@ -18,6 +18,7 @@ function oo_ = variance_decomposition_mc_analysis(NumberOfSimulations,type,dname % mh_conf_sig [double] 2 by 1 vector with upper % and lower bound of HPD intervals % oo_ [structure] Dynare structure where the results are saved. +% options_ [structure] Dynare options structure % % OUTPUTS % oo_ [structure] Dynare structure where the results are saved. @@ -93,7 +94,7 @@ if t3<1.0e-12 end if abs(t1-1)<1.0e-12 t1 = 1; - end + end p_mean = t1; p_median = t1; p_var = 0; @@ -101,8 +102,13 @@ if t3<1.0e-12 p_deciles = NaN(9,1); density = NaN; else - [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ... - posterior_moments(tmp,1,mh_conf_sig); + if options_.estimation.moments_posterior_density.indicator + [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ... + posterior_moments(tmp,1,mh_conf_sig); + else + [p_mean, p_median, p_var, hpd_interval, p_deciles] = ... + posterior_moments(tmp,0,mh_conf_sig); + end end eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.Mean.' name ' = p_mean;']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.Median.' name ' = p_median;']); @@ -110,4 +116,6 @@ eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.Variance.' name eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.HPDinf.' name ' = hpd_interval(1);']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.HPDsup.' name ' = hpd_interval(2);']); eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.deciles.' name ' = p_deciles;']); -eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.density.' name ' = density;']); \ No newline at end of file +if options_.estimation.moments_posterior_density.indicator + eval(['oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.density.' name ' = density;']); +end \ No newline at end of file diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index 1b99af33a..925e7c14a 100644 --- a/preprocessor/DynareBison.yy +++ b/preprocessor/DynareBison.yy @@ -114,7 +114,7 @@ class ParsingDriver; %token NOGRAPH NOMOMENTS NOPRINT NORMAL_PDF SAVE_DRAWS %token OBSERVATION_TRENDS OPTIM OPTIM_WEIGHTS ORDER OSR OSR_PARAMS MAX_DIM_COVA_GROUP ADVANCED OUTFILE OUTVARS OVERWRITE %token PARALLEL_LOCAL_FILES PARAMETERS PARAMETER_SET PARTIAL_INFORMATION PERFECT_FORESIGHT PERIODS PERIOD PLANNER_OBJECTIVE PLOT_CONDITIONAL_FORECAST PLOT_PRIORS PREFILTER PRESAMPLE -%token PERFECT_FORESIGHT_SETUP PERFECT_FORESIGHT_SOLVER POSTERIOR_KERNEL_DENSITY +%token PERFECT_FORESIGHT_SETUP PERFECT_FORESIGHT_SOLVER NO_POSTERIOR_KERNEL_DENSITY %token PRINT PRIOR_MC PRIOR_TRUNC PRIOR_MODE PRIOR_MEAN POSTERIOR_MODE POSTERIOR_MEAN POSTERIOR_MEDIAN PRUNING %token QUOTED_STRING %token QZ_CRITERIUM QZ_ZERO_THRESHOLD FULL DSGE_VAR DSGE_VARLAG DSGE_PRIOR_WEIGHT TRUNCATE @@ -1739,7 +1739,7 @@ estimation_options : o_datafile | o_silent_optimizer | o_proposal_distribution | o_student_degrees_of_freedom - | o_posterior_kernel_density + | o_no_posterior_kernel_density ; list_optim_option : QUOTED_STRING COMMA QUOTED_STRING @@ -2692,8 +2692,8 @@ o_mh_jscale : MH_JSCALE EQUAL non_negative_number { driver.option_num("mh_jscale o_optim : OPTIM EQUAL '(' optim_options ')'; o_tarb_optim : TARB_OPTIM EQUAL '(' tarb_optim_options ')'; o_proposal_distribution : PROPOSAL_DISTRIBUTION EQUAL symbol { driver.option_str("proposal_distribution", $3); }; -o_posterior_kernel_density : POSTERIOR_KERNEL_DENSITY - { driver.option_num("posterior_kernel_density.indicator", "1"); } +o_no_posterior_kernel_density : NO_POSTERIOR_KERNEL_DENSITY + { driver.option_num("moments_posterior_density.indicator", "0"); } ; o_student_degrees_of_freedom : STUDENT_DEGREES_OF_FREEDOM EQUAL INT_NUMBER { driver.option_num("student_degrees_of_freedom", $3); }; o_mh_init_scale : MH_INIT_SCALE EQUAL non_negative_number { driver.option_num("mh_init_scale", $3); }; diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 1445c44db..82209584a 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -377,7 +377,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2 montecarlo {return token::MONTECARLO;} distribution_approximation {return token::DISTRIBUTION_APPROXIMATION;} proposal_distribution {return token::PROPOSAL_DISTRIBUTION;} -posterior_kernel_density {return token::POSTERIOR_KERNEL_DENSITY;} +no_posterior_kernel_density {return token::NO_POSTERIOR_KERNEL_DENSITY;} student_degrees_of_freedom {return token::STUDENT_DEGREES_OF_FREEDOM;} alpha { diff --git a/tests/estimation/fs2000.mod b/tests/estimation/fs2000.mod index 3b51be014..2eaeff302 100644 --- a/tests/estimation/fs2000.mod +++ b/tests/estimation/fs2000.mod @@ -82,4 +82,4 @@ varobs gp_obs gy_obs; options_.solve_tolf = 1e-12; -estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=3000,mh_nblocks=2,mh_jscale=0.8,moments_varendo,selected_variables_only) y m; +estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=3000,mh_nblocks=2,mh_jscale=0.8,moments_varendo,selected_variables_only,contemporaneous_correlation,smoother,forecast=8) y m; From 7ebbc41d8af472342f813b9891c1a5c992e96c44 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 12 Oct 2015 15:20:15 +0200 Subject: [PATCH 3/3] Get rid of eval statements by using indirect indexing of structures --- matlab/GetPosteriorParametersStatistics.m | 2 +- matlab/PosteriorIRF.m | 24 ++++----- ...ional_variance_decomposition_mc_analysis.m | 22 ++++---- matlab/correlation_mc_analysis.m | 42 +++++++-------- matlab/covariance_mc_analysis.m | 54 +++++++++---------- matlab/variance_decomposition_mc_analysis.m | 20 +++---- 6 files changed, 81 insertions(+), 83 deletions(-) diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m index f9bfedbd0..c51e0a2ba 100644 --- a/matlab/GetPosteriorParametersStatistics.m +++ b/matlab/GetPosteriorParametersStatistics.m @@ -389,6 +389,6 @@ oo.posterior_density.(type).(name) = density; function [post_mean,hpd_interval,post_var] = Extractoo(oo,name,type) hpd_interval = zeros(2,1); post_mean = oo.posterior_mean.(type).(name); -hpd_interval(1) = oo.posterior_hpdinf.(type).(name); +hpd_interval(1) = oo.posterior_hpdinf.(type).(name); hpd_interval(2) = oo.posterior_hpdsup.(type).(name); post_var = oo.posterior_variance.(type).(name); \ No newline at end of file diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m index 8c094a24f..be4246100 100644 --- a/matlab/PosteriorIRF.m +++ b/matlab/PosteriorIRF.m @@ -304,12 +304,12 @@ clear STOCK_IRF_DSGE; for i = irf_shocks_indx for j = 1:nvar name = [deblank(M_.endo_names(IndxVariables(j),:)) '_' deblank(tit(i,:))]; - eval(['oo_.PosteriorIRF.dsge.Mean.' name ' = MeanIRF(:,j,i);']); - eval(['oo_.PosteriorIRF.dsge.Median.' name ' = MedianIRF(:,j,i);']); - eval(['oo_.PosteriorIRF.dsge.Var.' name ' = VarIRF(:,j,i);']); - eval(['oo_.PosteriorIRF.dsge.deciles.' name ' = DistribIRF(:,:,j,i);']); - eval(['oo_.PosteriorIRF.dsge.HPDinf.' name ' = HPDIRF(:,1,j,i);']); - eval(['oo_.PosteriorIRF.dsge.HPDsup.' name ' = HPDIRF(:,2,j,i);']); + oo_.PosteriorIRF.dsge.Mean.(name) = MeanIRF(:,j,i); + oo_.PosteriorIRF.dsge.Median.(name) = MedianIRF(:,j,i); + oo_.PosteriorIRF.dsge.Var.(name) = VarIRF(:,j,i); + oo_.PosteriorIRF.dsge.deciles.(name) = DistribIRF(:,:,j,i); + oo_.PosteriorIRF.dsge.HPDinf.(name) = HPDIRF(:,1,j,i); + oo_.PosteriorIRF.dsge.HPDsup.(name) = HPDIRF(:,2,j,i); end end @@ -341,12 +341,12 @@ if MAX_nirfs_dsgevar for i = irf_shocks_indx for j = 1:nvar name = [deblank(M_.endo_names(IndxVariables(j),:)) '_' deblank(tit(i,:))]; - eval(['oo_.PosteriorIRF.bvardsge.Mean.' name ' = MeanIRFdsgevar(:,j,i);']); - eval(['oo_.PosteriorIRF.bvardsge.Median.' name ' = MedianIRFdsgevar(:,j,i);']); - eval(['oo_.PosteriorIRF.bvardsge.Var.' name ' = VarIRFdsgevar(:,j,i);']); - eval(['oo_.PosteriorIRF.bvardsge.deciles.' name ' = DistribIRFdsgevar(:,:,j,i);']); - eval(['oo_.PosteriorIRF.bvardsge.HPDinf.' name ' = HPDIRFdsgevar(:,1,j,i);']); - eval(['oo_.PosteriorIRF.bvardsge.HPDsup.' name ' = HPDIRFdsgevar(:,2,j,i);']); + oo_.PosteriorIRF.bvardsge.Mean.(name) = MeanIRFdsgevar(:,j,i); + oo_.PosteriorIRF.bvardsge.Median.(name) = MedianIRFdsgevar(:,j,i); + oo_.PosteriorIRF.bvardsge.Var.(name) = VarIRFdsgevar(:,j,i); + oo_.PosteriorIRF.bvardsge.deciles.(name) = DistribIRFdsgevar(:,:,j,i); + oo_.PosteriorIRF.bvardsge.HPDinf.(name) = HPDIRFdsgevar(:,1,j,i); + oo_.PosteriorIRF.bvardsge.HPDsup.(name) = HPDIRFdsgevar(:,2,j,i); end end end diff --git a/matlab/conditional_variance_decomposition_mc_analysis.m b/matlab/conditional_variance_decomposition_mc_analysis.m index 2c724083a..c165a62be 100644 --- a/matlab/conditional_variance_decomposition_mc_analysis.m +++ b/matlab/conditional_variance_decomposition_mc_analysis.m @@ -62,11 +62,11 @@ end name = [ var_list(endogenous_variable_index,:) '.' exo ]; if isfield(oo_, [ TYPE 'TheoreticalMoments' ]) - eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments;']) + temporary_structure = oo_.([TYPE 'TheoreticalMoments']); if isfield(temporary_structure,'dsge') - eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge;']) + temporary_structure = oo_.([TYPE 'TheoreticalMoments']).dsge; if isfield(temporary_structure,'ConditionalVarianceDecomposition') - eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.Mean;']) + temporary_structure = oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.Mean; if isfield(temporary_structure,name) if sum(Steps-temporary_structure.(name)(1,:)) == 0 % Nothing (new) to do here... @@ -112,13 +112,13 @@ for i=1:length(Steps) p_hpdinf(i) = hpd_interval(1); p_hpdsup(i) = hpd_interval(2); end -eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.Steps = Steps;']); -eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.Mean.' name ' = p_mean;']); -eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.Median.' name ' = p_median;']); -eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.Variance.' name ' = p_variance;']); -eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.HPDinf.' name ' = p_hpdinf;']); -eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.HPDsup.' name ' = p_hpdsup;']); -eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.deciles.' name ' = p_deciles;']); +oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.Steps = Steps; +oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.Mean.(name) = p_mean; +oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.Median.(name) = p_median; +oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.Variance.(name) = p_variance; +oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.HPDinf.(name) = p_hpdinf; +oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.HPDsup.(name) = p_hpdsup; +oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.deciles.(name) = p_deciles; if options_.estimation.moments_posterior_density.indicator - eval(['oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.density.' name ' = p_density;']); + oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecomposition.density.(name) = p_density; end \ No newline at end of file diff --git a/matlab/correlation_mc_analysis.m b/matlab/correlation_mc_analysis.m index 28a5c1af7..def6c94a9 100644 --- a/matlab/correlation_mc_analysis.m +++ b/matlab/correlation_mc_analysis.m @@ -44,15 +44,15 @@ else end if isfield(oo_,[TYPE 'TheoreticalMoments']) - eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments;']) + temporary_structure = oo_.([TYPE, 'TheoreticalMoments']); if isfield(temporary_structure,'dsge') - eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge;']) + temporary_structure = oo_.([TYPE, 'TheoreticalMoments']).dsge; if isfield(temporary_structure,'correlation') - eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge.correlation.Mean;']) + temporary_structure = oo_.([TYPE, 'TheoreticalMoments']).dsge.correlation.Mean; if isfield(temporary_structure,deblank(var1)) - eval(['temporary_structure_1 = oo_.' TYPE 'TheoreticalMoments.dsge.correlation.Mean.' var1 ';']) + temporary_structure_1 = oo_.([TYPE, 'TheoreticalMoments']).dsge.correlation.Mean.(var1); if isfield(temporary_structure_1,deblank(var2)) - eval(['temporary_structure_2 = temporary_structure_1.' var2 ';']) + temporary_structure_2 = temporary_structure_1.(var2); l1 = length(temporary_structure_2); if l1