From b3877b3a9e86a2441883c3f0b95fbf1284bfa4d5 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Tue, 14 Jun 2016 15:09:05 +0200 Subject: [PATCH 1/2] Save chain's proposal density and make load_mh_file and mh_recover load it Closes #504 --- matlab/posterior_sampler.m | 2 +- matlab/posterior_sampler_core.m | 3 +++ matlab/posterior_sampler_initialization.m | 27 ++++++++++++++++++++--- tests/estimation/fs2000.mod | 19 +++++++++++++++- 4 files changed, 46 insertions(+), 5 deletions(-) diff --git a/matlab/posterior_sampler.m b/matlab/posterior_sampler.m index 9efabb4d8..3e187a19e 100644 --- a/matlab/posterior_sampler.m +++ b/matlab/posterior_sampler.m @@ -61,7 +61,7 @@ objective_function_penalty_base = Inf; vv = sampler_options.invhess; % Initialization of the sampler -[ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ... +[ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d, bayestopt_] = ... posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2); diff --git a/matlab/posterior_sampler_core.m b/matlab/posterior_sampler_core.m index ec749de63..c677e02cc 100644 --- a/matlab/posterior_sampler_core.m +++ b/matlab/posterior_sampler_core.m @@ -110,6 +110,9 @@ end sampler_options.xparam1 = xparam1; if ~isempty(d), sampler_options.proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale); + %store information for load_mh_file + record.ProposalCovariance=d; + record.ProposalScaleVec=bayestopt_.jscale; end block_iter=0; diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m index 75df703c8..3adca10a5 100644 --- a/matlab/posterior_sampler_initialization.m +++ b/matlab/posterior_sampler_initialization.m @@ -1,6 +1,6 @@ -function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d ] = ... +function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d, bayestopt_] = ... posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) -% function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d ] = ... +% function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d, bayestopt_] = ... % metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) % Metropolis-Hastings initialization. % @@ -33,6 +33,7 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npa % draws % o MAX_nruns [scalar] maximum number of draws in each MH-file on the harddisk % o d [double] (p*p) matrix, Cholesky decomposition of the posterior covariance matrix (at the mode). +% o bayestopt_ [structure] estimation options structure % % SPECIAL REQUIREMENTS % None. @@ -224,6 +225,7 @@ if ~options_.load_mh_file && ~options_.mh_recover record.LastFileNumber = AnticipatedNumberOfFiles ; record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile; record.MCMCConcludedSuccessfully = 0; + record.MCMC_sampler=options_.posterior_sampler_options.posterior_sampling_method; fprintf('Ok!\n'); id = write_mh_history_file(MetropolisFolder, ModelName, record); disp(['Estimation::mcmc: Details about the MCMC are available in ' BaseName '_mh_history_' num2str(id) '.mat']) @@ -283,6 +285,7 @@ elseif options_.load_mh_file && ~options_.mh_recover end ilogpo2 = record.LastLogPost; ix2 = record.LastParameters; + [d,bayestopt_]=set_proposal_density_to_previous_value(record,options_,bayestopt_); FirstBlock = 1; NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1); fprintf('Estimation::mcmc: I am writing a new mh-history file... '); @@ -363,7 +366,7 @@ elseif options_.mh_recover FirstLine = ones(NumberOfBlocks,1); LastFileFullIndicator=1; end - + [d,bayestopt_]=set_proposal_density_to_previous_value(record,options_,bayestopt_); %% Now find out what exactly needs to be redone % 1. Check if really something needs to be done % How many mh files should we have ? @@ -447,3 +450,21 @@ elseif options_.mh_recover write_mh_history_file(MetropolisFolder, ModelName, record); end end + +function [d,bayestopt_]=set_proposal_density_to_previous_value(record,options_,bayestopt_) + if isfield(record,'ProposalCovariance') && isfield(record,'ProposalCovariance') + if isfield(record,'MCMC_sampler') + if ~strcmp(record.MCMC_sampler,options_.posterior_sampler_options.posterior_sampling_method) + error(fprintf('Estimation::mcmc: The posterior_sampling_method differs from the one of the original chain. Please reset it to %s',record.MCMC_sampler)) + end + end + fprintf('Estimation::mcmc: Recovering the previous proposal density\n.') + d=record.ProposalCovariance; + bayestopt_.jscale=record.ProposalScaleVec; + else + if options_.mode_compute~=0 + fprintf('Estimation::mcmc: No stored previous proposal density found, continuing with the one implied by mode_compute\n.'); + elseif ~isempty(options_.mode_file) + fprintf('Estimation::mcmc: No stored previous proposal density found, continuing with the one implied by the mode_file\n.'); + end + end diff --git a/tests/estimation/fs2000.mod b/tests/estimation/fs2000.mod index 2eaeff302..640f1961f 100644 --- a/tests/estimation/fs2000.mod +++ b/tests/estimation/fs2000.mod @@ -82,4 +82,21 @@ 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,contemporaneous_correlation,smoother,forecast=8) y m; +estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=3000,mh_nblocks=1,mh_jscale=0.8,moments_varendo,selected_variables_only,contemporaneous_correlation,smoother,forecast=8) y m; +%test load_mh_file option +options_.smoother=0; +options_.moments_varendo=0; +options_.forecast=0; +copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat'],[M_.dname '_mh1_blck1.mat']) +estimation(mode_compute=0,mode_file=fs2000_mode,order=1, datafile=fsdat_simul, nobs=192, loglinear, mh_replic=1500, mh_nblocks=1, mh_jscale=0.8); +hh=eye(size(bayestopt_.name,1)); +save('fs2000_mode','hh','-append') +estimation(mode_compute=0,mode_file=fs2000_mode,order=1, datafile=fsdat_simul, nobs=192, loglinear, mh_replic=1500, mh_nblocks=1, mh_jscale=10,load_mh_file); + + +temp1=load([M_.dname '_mh1_blck1.mat']); +temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat']); + +if max(max(abs(temp1.x2-temp2.x2)))>1e-10 + error('Draws of unaffected chain are not the same') +end From 1084a2d452a7f218c5ef43551a71957a56359c08 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Tue, 14 Jun 2016 15:18:05 +0200 Subject: [PATCH 2/2] Document new behavior of mh_recover and mh_replic --- doc/dynare.texi | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index 6a938c3df..a2ace0642 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -5124,8 +5124,9 @@ Metropolis-Hastings chain. Default: 2*@code{mh_scale} @anchor{mh_recover} Attempts to recover a Metropolis-Hastings simulation that crashed prematurely, starting with the last available saved @code{mh}-file. Shouldn't be used together with -@code{load_mh_file} or a different @code{mh_replic} than in the crashed run. -To assure a neat continuation of the chain with the same proposal density, you should +@code{load_mh_file} or a different @code{mh_replic} than in the crashed run. Since Dynare 4.5 +the proposal density from the previous run will automatically be loaded. In older versions, +to assure a neat continuation of the chain with the same proposal density, you should provide the @code{mode_file} used in the previous run or the same user-defined @code{mcmc_jumping_covariance} when using this option. @@ -5293,7 +5294,12 @@ when finite values are required for computational reasons. Default: @code{1e7} @item load_mh_file @anchor{load_mh_file} Tells Dynare to add to previous Metropolis-Hastings simulations instead of starting from -scratch. Shouldn't be used together with @code{mh_recover} +scratch. Since Dynare 4.5 +the proposal density from the previous run will automatically be loaded. In older versions, +to assure a neat continuation of the chain with the same proposal density, you should +provide the @code{mode_file} used in the previous +run or the same user-defined @code{mcmc_jumping_covariance} when using this option. +Shouldn't be used together with @code{mh_recover}. @item optim = (@var{NAME}, @var{VALUE}, ...) @anchor{optim}