From 63986a0ebf6478212e6d5ba70bd51f47f13510b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Scylla=29?= Date: Wed, 18 Dec 2013 16:39:41 +0100 Subject: [PATCH] Closes #567. --- doc/dynare.texi | 18 +++++++++++++----- matlab/PosteriorIRF.m | 14 +++++--------- matlab/dynare_estimation_1.m | 11 ++++++++++- matlab/global_initialization.m | 2 +- matlab/metropolis_draw.m | 16 +++++++++++++++- matlab/posterior_analysis.m | 2 +- matlab/prior_posterior_statistics.m | 19 ++----------------- preprocessor/DynareBison.yy | 4 +++- preprocessor/DynareFlex.ll | 1 + 9 files changed, 51 insertions(+), 36 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index ba6d7bcc8..4a312ee6e 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4494,12 +4494,20 @@ algorithm. For the time being, @code{mh_replic} should be larger than @code{1200}. Default: @code{20000} @item sub_draws = @var{INTEGER} -@anchor{sub_draws} number of draws from the Metropolis iterations that -are used to compute posterior distribution of various objects (smoothed -variable, smoothed shocks, forecast, moments, IRF). @code{sub_draws} should be smaller than -the total number of Metropolis draws available. Default: -@code{min(1200,0.25*Total number of draws)} +@anchor{sub_draws} number of draws from the MCMC that are used to +compute posterior distribution of various objects (smoothed variable, +smoothed shocks, forecast, moments, IRF). The draws used to compute +these posterior moments are sampled uniformly in the estimated empirical +posterior distribution (@i{ie} draws of the MCMC). @code{sub_draws} +should be smaller than the total number of MCMC draws available. +Default: @code{min(posterior_max_subsample_draws,0.25*Total number of +draws)} +@item posterior_max_subsample_draws = @var{INTEGER} +@anchor{posterior_max_subsample_draws} maximum number of draws from the +MCMC used to compute posterior distribution of various objects (smoothed +variable, smoothed shocks, forecast, moments, IRF), if not overriden by +option @ref{sub_draws}. @item mh_nblocks = @var{INTEGER} @anchor{mh_nblocks} Number of parallel chains for Metropolis-Hastings algorithm. Default: diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m index 018c6ad3a..cf9737f0f 100644 --- a/matlab/PosteriorIRF.m +++ b/matlab/PosteriorIRF.m @@ -95,9 +95,8 @@ end delete([MhDirectoryName filesep M_.fname '_IRF_DSGEs*.mat']); delete([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs*.mat']); if strcmpi(type,'posterior') - load_last_mh_history_file(MhDirectoryName, M_.fname); - TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); - NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); + B = options_.sub_draws; + options_.B = B; elseif strcmpi(type,'gsa') RootDirectoryName = CheckPath('gsa',M_.dname); if options_.opt_gsa.pprior @@ -107,13 +106,10 @@ elseif strcmpi(type,'gsa') end x=[lpmat0(istable,:) lpmat(istable,:)]; clear lpmat istable - NumberOfDraws=size(x,1); - B=NumberOfDraws; options_.B = B; + B=size(x,1); options_.B = B; else% type = 'prior' - NumberOfDraws = options_.prior_draws; -end -if ~strcmpi(type,'gsa') - B = min([round(.5*NumberOfDraws),500]); options_.B = B; + B = options_.prior_draws; + options_.B = B; end try delete([MhDirectoryName filesep M_.fname '_irf_dsge*.mat']) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 0578791d7..d925acc58 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -827,14 +827,23 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... else load([M_.fname '_results'],'oo_'); end - metropolis_draw(1); + error_flag = metropolis_draw(1); if options_.bayesian_irf + if error_flag + error('Estimation::mcmc: I cannot compute the posterior IRFs!') + end PosteriorIRF('posterior'); end if options_.moments_varendo + if error_flag + error('Estimation::mcmc: I cannot compute the posterior moments for the endogenous variables!') + end oo_ = compute_moments_varendo('posterior',options_,M_,oo_,var_list_); end if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast + if error_flag + error('Estimation::mcmc: I cannot compute the posterior statistics!') + end prior_posterior_statistics('posterior',dataset_); end xparam = get_posterior_parameters('mean'); diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index de4741c13..04b2242a2 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -355,7 +355,6 @@ options_.dataset.xls_range = NaN; options_.Harvey_scale_factor = 10; options_.MaxNumberOfBytes = 1e6; options_.MaximumNumberOfMegaBytes = 111; -options_.PosteriorSampleSize = 1000; options_.analytic_derivation = 0; options_.analytic_derivation_mode = 0; options_.bayesian_irf = 0; @@ -399,6 +398,7 @@ options_.presample = 0; options_.prior_trunc = 1e-10; options_.smoother = 0; options_.student_degrees_of_freedom = 3; +options_.posterior_max_subsample_draws = 1200; options_.sub_draws = []; options_.use_mh_covariance_matrix = 0; options_.gradient_method = 2; diff --git a/matlab/metropolis_draw.m b/matlab/metropolis_draw.m index 8e1f10453..771ba0a43 100644 --- a/matlab/metropolis_draw.m +++ b/matlab/metropolis_draw.m @@ -32,6 +32,9 @@ function [xparams, logpost]=metropolis_draw(init) global options_ estim_params_ M_ persistent mh_nblck NumberOfDraws BaseName FirstLine FirstMhFile MAX_nruns +xparams = 0; +logpost = 0; + if init nvx = estim_params_.nvx; nvn = estim_params_.nvn; @@ -39,7 +42,6 @@ if init ncn = estim_params_.ncn; np = estim_params_.np ; npar = nvx+nvn+ncx+ncn+np; - MetropolisFolder = CheckPath('metropolis',M_.dname); FileName = M_.fname; BaseName = [MetropolisFolder filesep FileName]; @@ -52,6 +54,18 @@ if init NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); mh_nblck = options_.mh_nblck; + % set sub_draws option if empty + if isempty(options_.sub_draws) + options_.sub_draws = min(options_.posterior_max_subsample_draws, round(.25*NumberOfDraws)); + else + if options_.sub_draws>NumberOfDraws + skipline() + disp(['Estimation::mcmc: The value of option sub_draws (' num2str(options_.sub_draws) ') is greater than the number of available draws in the MCMC (' num2str(NumberOfDraws) ')!']) + disp('Estimation::mcmc: You can either change the value of sub_draws, reduce the value of mh_drop, or run another mcmc (with the load_mh_file option).') + skipline() + xparams = 1; % xparams is interpreted as an error flag + end + end return end diff --git a/matlab/posterior_analysis.m b/matlab/posterior_analysis.m index df978c647..0382f923f 100644 --- a/matlab/posterior_analysis.m +++ b/matlab/posterior_analysis.m @@ -17,7 +17,7 @@ function oo_ = posterior_analysis(type,arg1,arg2,arg3,options_,M_,oo_) % along with Dynare. If not, see . info = check_posterior_analysis_data(type,M_); -SampleSize = options_.PosteriorSampleSize; +SampleSize = options_.sub_draws; switch info case 0 disp('check_posterior_analysis_data:: Can''t find any mcmc file!') diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m index fa4ff29ad..f9192164e 100644 --- a/matlab/prior_posterior_statistics.m +++ b/matlab/prior_posterior_statistics.m @@ -69,21 +69,7 @@ maxlag = M_.maximum_endo_lag; if strcmpi(type,'posterior') DirectoryName = CheckPath('metropolis',M_.dname); - load_last_mh_history_file(DirectoryName, M_.fname); - FirstMhFile = record.KeepedDraws.FirstMhFile; - FirstLine = record.KeepedDraws.FirstLine; - TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); LastMhFile = TotalNumberOfMhFiles; - TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); - NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); - clear record; - if ~isempty(options_.sub_draws) - B = options_.sub_draws; - if B > NumberOfDraws - B = NumberOfDraws; - end - else - B = min(1200, round(0.25*NumberOfDraws)); - end + B = options_.sub_draws; elseif strcmpi(type,'gsa') RootDirectoryName = CheckPath('gsa',M_.dname); if options_.opt_gsa.pprior @@ -99,8 +85,7 @@ elseif strcmpi(type,'gsa') x=lpmat(istable,:); end clear lpmat lpmat0 istable - NumberOfDraws=size(x,1); - B=NumberOfDraws; + B = size(x,1); elseif strcmpi(type,'prior') DirectoryName = CheckPath('prior',M_.dname); B = options_.prior_draws; diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index ea7d82e0d..3e95184c9 100644 --- a/preprocessor/DynareBison.yy +++ b/preprocessor/DynareBison.yy @@ -110,7 +110,7 @@ class ParsingDriver; %token KALMAN_ALGO KALMAN_TOL SUBSAMPLES OPTIONS TOLF %token LAPLACE LIK_ALGO LIK_INIT LINEAR LOAD_IDENT_FILES LOAD_MH_FILE LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR LYAPUNOV %token LYAPUNOV_FIXED_POINT_TOL LYAPUNOV_DOUBLING_TOL LYAPUNOV_SQUARE_ROOT_SOLVER_TOL LOG_DEFLATOR LOG_TREND_VAR LOG_GROWTH_FACTOR MARKOWITZ MARGINAL_DENSITY MAX MAXIT -%token MFS MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER MIN MINIMAL_SOLVING_PERIODS +%token MFS MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER POSTERIOR_MAX_SUBSAMPLE_DRAWS MIN MINIMAL_SOLVING_PERIODS %token MODE_CHECK MODE_CHECK_NEIGHBOURHOOD_SIZE MODE_CHECK_SYMMETRIC_PLOTS MODE_CHECK_NUMBER_OF_POINTS MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS ABS SIGN %token MODEL_DIAGNOSTICS MODIFIEDHARMONICMEAN MOMENTS_VARENDO DIFFUSE_FILTER SUB_DRAWS TAPER_STEPS GEWEKE_INTERVAL MCMC_JUMPING_COVARIANCE %token NAME @@ -1633,6 +1633,7 @@ estimation_options : o_datafile | o_geweke_interval | o_mcmc_jumping_covariance | o_irf_plot_threshold + | o_posterior_max_subsample_draws ; list_optim_option : QUOTED_STRING COMMA QUOTED_STRING @@ -2423,6 +2424,7 @@ o_subsample_name : symbol EQUAL date_expr ':' date_expr ; o_conf_sig : CONF_SIG EQUAL non_negative_number { driver.option_num("conf_sig", $3); }; o_mh_replic : MH_REPLIC EQUAL INT_NUMBER { driver.option_num("mh_replic", $3); }; +o_posterior_max_subsample_draws : POSTERIOR_MAX_SUBSAMPLE_DRAWS EQUAL INT_NUMBER { driver.option_num("posterior_max_subsample_draws", $3); }; o_mh_drop : MH_DROP EQUAL non_negative_number { driver.option_num("mh_drop", $3); }; o_mh_jscale : MH_JSCALE EQUAL non_negative_number { driver.option_num("mh_jscale", $3); }; o_optim : OPTIM EQUAL '(' optim_options ')'; diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 7343a5ad1..9b6b5ed46 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -284,6 +284,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 dsge_var {return token::DSGE_VAR;} dsge_varlag {return token::DSGE_VARLAG;} moments_varendo {return token::MOMENTS_VARENDO;} +posterior_max_subsample_draws {return token::POSTERIOR_MAX_SUBSAMPLE_DRAWS;} filtered_vars {return token::FILTERED_VARS;} filter_step_ahead {return token::FILTER_STEP_AHEAD;} relative_irf {return token::RELATIVE_IRF;}