stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'rattoma/dynare-mh_initialize_from_previous_mcmc' 

Ref. !1841
time-shift
Sébastien Villemot 2021-05-11 18:51:57 +02:00
parent 306e963efa a0f2b2399b
commit 0d6bc47158
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
5 changed files with 222 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
doc/manual/source/the-model-file.rst
View File

@ -5557,6 +5557,48 @@ block decomposition of the model (see :opt:`block`).
density, and posterior statistics from an existing ``_results``
file instead of recomputing them.
.. option:: mh_initialize_from_previous_mcmc
This option allows to pick initial values for new MCMC from a previous one,
where the model specification, the number of estimated parameters,
(some) prior might have changed (so a situation where ``load_mh_file`` would not work).
If an additional parameter is estimated, it is automatically initialized from prior_draw.
Note that, if this option is used to skip the optimization step, you should use a sampling method which does not require
a proposal density, like slice. Otherwise, optimization should always be done beforehand or a mode file with
an appropriate posterior covariance matrix should be used.
.. option:: mh_initialize_from_previous_mcmc_directory = FILENAME
If ``mh_initialize_from_previous_mcmc`` is set, users must provide here
the path to the standard FNAME folder from where to load prior definitions and
last MCMC values to be used to initialize the new MCMC.
Example: if previous project directory is ``/my_previous_dir`` and FNAME is ``mymodel``,
users should set the option as
``mh_initialize_from_previous_mcmc_directory = '/my_previous_dir/mymodel'``
Dynare will then look for the last record file into
``/my_previous_dir/mymodel/metropolis/mymodel_mh_history_<LAST>.mat``
and for the prior definition file into
``/my_previous_dir/mymodel/prior/definition.mat``
.. option:: mh_initialize_from_previous_mcmc_record = FILENAME
If ``mh_initialize_from_previous_mcmc`` is set, and whenever the standard file or directory tree
is not applicable to load initial values, users may directly provide here
the path to the record file from which to load
values to be used to initialize the new MCMC.
.. option:: mh_initialize_from_previous_mcmc_prior = FILENAME
If ``mh_initialize_from_previous_mcmc`` is set, and whenever the standard file or directory tree
is not applicable to load initial values, users may directly provide here
the path to the prior definition file, to get info in the priors used in previous MCMC.
.. option:: optim = (NAME, VALUE, ...)
A list of NAME and VALUE pairs. Can be used to set options for

4
matlab/default_option_values.m
View File

@ -396,6 +396,10 @@ options_.mh_tune_jscale.c1 = .02;
options_.mh_tune_jscale.c2 = .05;
options_.mh_tune_jscale.c3 = 4;
options_.mh_init_scale = 2*options_.mh_jscale;
options_.mh_initialize_from_previous_mcmc.status = false;
options_.mh_initialize_from_previous_mcmc.directory = '';
options_.mh_initialize_from_previous_mcmc.record = '';
options_.mh_initialize_from_previous_mcmc.prior = '';
options_.mh_mode = 1;
options_.mh_nblck = 2;
options_.mh_recover = false;

90
matlab/posterior_sampler_initialization.m
View File

@ -115,13 +115,80 @@ if ~options_.load_mh_file && ~options_.mh_recover
if isempty(d)
prior_draw(bayestopt_,options_.prior_trunc);
end
if options_.mh_initialize_from_previous_mcmc.status
PreviousFolder0 = options_.mh_initialize_from_previous_mcmc.directory;
RecordFile0 = options_.mh_initialize_from_previous_mcmc.record;
PriorFile0 = options_.mh_initialize_from_previous_mcmc.prior;
if ~isempty(RecordFile0)
%% check for proper filesep char in user defined paths
RecordFile0=strrep(RecordFile0,'\',filesep);
if isempty(dir(RecordFile0))
disp('Estimation::mcmc: wrong value for mh_initialize_from_previous_mcmc_record option')
error('Estimation::mcmc: path to record file is not found')
else
record0=load(RecordFile0);
end
record0=record0.record;
MetropolisFolder0 = fileparts(RecordFile0);
PreviousFolder0=fileparts(MetropolisFolder0);
[~, ModelName0]=fileparts(PreviousFolder0);
else
%% check for proper filesep char in user defined paths
PreviousFolder0=strrep(PreviousFolder0,'\',filesep);
MetropolisFolder0 = [PreviousFolder0 filesep 'metropolis'];
[~, ModelName0]=fileparts(PreviousFolder0);
record0=load_last_mh_history_file(MetropolisFolder0, ModelName0);
end
if ~isnan(record0.MCMCConcludedSuccessfully) && ~record0.MCMCConcludedSuccessfully
error('Estimation::mcmc: You are trying to load an MCMC that did not finish successfully. Please use mh_recover.')
end
% mh_files = dir([ MetropolisFolder0 filesep ModelName0 '_mh*.mat']);
% if ~length(mh_files)
% error('Estimation::mcmc: I cannot find any MH file to load here!')
% end
disp('Estimation::mcmc: Initializing from past Metropolis-Hastings simulations...')
disp(['Estimation::mcmc: Past MH path ' MetropolisFolder0 ])
disp(['Estimation::mcmc: Past model name ' ModelName0 ])
fprintf(fidlog,' Loading initial values from previous MH\n');
fprintf(fidlog,' Past MH path: %s\n', MetropolisFolder0 );
fprintf(fidlog,' Past model name: %s\n', ModelName0);
fprintf(fidlog,' \n');
past_number_of_blocks = record0.Nblck;
if past_number_of_blocks ~= NumberOfBlocks
disp('Estimation::mcmc: The specified number of blocks doesn''t match with the previous number of blocks!')
disp(['Estimation::mcmc: You declared ' int2str(NumberOfBlocks) ' blocks, but the previous number of blocks was ' int2str(past_number_of_blocks) '.'])
disp(['Estimation::mcmc: I will run the Metropolis-Hastings with ' int2str(past_number_of_blocks) ' blocks.' ])
NumberOfBlocks = past_number_of_blocks;
options_.mh_nblck = NumberOfBlocks;
end
if ~isempty(PriorFile0)
if isempty(dir(PriorFile0))
disp('Estimation::mcmc: wrong value for mh_initialize_from_previous_mcmc_prior option')
error('Estimation::mcmc: path to prior file is not found')
else
bayestopt0 = load(PriorFile0);
end
else
bayestopt0 = load([PreviousFolder0 filesep 'prior' filesep 'definition.mat']);
end
[common_parameters,IA,IB] = intersect(bayestopt_.name,bayestopt0.bayestopt_.name);
new_estimated_parameters = ~ismember(bayestopt_.name,bayestopt0.bayestopt_.name);
ix2 = zeros(NumberOfBlocks,npar);
ilogpo2 = zeros(NumberOfBlocks,1);
ilogpo2(:,1) = record0.LastLogPost;
ix2(:,IA) = record0.LastParameters(:,IB);
else
new_estimated_parameters = true(1,npar);
end
% Find initial values for the NumberOfBlocks chains...
if NumberOfBlocks > 1% Case 1: multiple chains
if NumberOfBlocks > 1 || options_.mh_initialize_from_previous_mcmc.status% Case 1: multiple chains
set_dynare_seed('default');
fprintf(fidlog,[' Initial values of the parameters:\n']);
disp('Estimation::mcmc: Searching for initial values...')
ix2 = zeros(NumberOfBlocks,npar);
ilogpo2 = zeros(NumberOfBlocks,1);
if ~options_.mh_initialize_from_previous_mcmc.status
ix2 = zeros(NumberOfBlocks,npar);
ilogpo2 = zeros(NumberOfBlocks,1);
end
for j=1:NumberOfBlocks
validate = 0;
init_iter = 0;
@ -133,7 +200,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
candidate = rand_multivariate_normal( transpose(xparam1), d * options_.mh_init_scale, npar);
end
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
ix2(j,:) = candidate;
ix2(j,new_estimated_parameters) = candidate(new_estimated_parameters);
ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
if ~isfinite(ilogpo2(j)) % if returned log-density is
% Inf or Nan (penalized value)
@ -212,11 +279,16 @@ if ~options_.load_mh_file && ~options_.mh_recover
record.MAX_nruns=MAX_nruns;
record.AcceptanceRatio = zeros(1,NumberOfBlocks);
record.FunctionEvalPerIteration = NaN(1,NumberOfBlocks);
for j=1:NumberOfBlocks
% we set a different seed for the random generator for each block then we record the corresponding random generator state (vector)
set_dynare_seed(options_.DynareRandomStreams.seed+j);
% record.Seeds keeps a vector of the random generator state and not the scalar seed despite its name
[record.InitialSeeds(j).Unifor,record.InitialSeeds(j).Normal] = get_dynare_random_generator_state();
if options_.mh_initialize_from_previous_mcmc.status
record.InitialSeeds = record0.LastSeeds;
else
for j=1:NumberOfBlocks
% we set a different seed for the random generator for each block then we record the corresponding random generator state (vector)
set_dynare_seed(options_.DynareRandomStreams.seed+j);
% record.Seeds keeps a vector of the random generator state and not the scalar seed despite its name
[record.InitialSeeds(j).Unifor,record.InitialSeeds(j).Normal] = get_dynare_random_generator_state();
end
end
record.InitialParameters = ix2;
record.InitialLogPost = ilogpo2;

3
tests/Makefile.am
View File

@ -44,6 +44,7 @@ MODFILES = \
estimation/fs2000_calibrated_covariance.mod \
estimation/fs2000_model_comparison.mod \
estimation/fs2000_fast.mod \
estimation/fs2000_init_from_previous.mod \
estimation/ls2003_endog_prior_restrict_estimation.mod \
estimation/independent_mh/fs2000_independent_mh.mod \
estimation/MH_recover/fs2000_recover.mod \
@ -647,6 +648,8 @@ AIM/ls2003_2L2L_AIM.o.trs: AIM/ls2003_2L2L.o.trs
estimation/fs2000_model_comparison.m.trs: estimation/fs2000.m.trs estimation/fs2000_initialize_from_calib.m.trs estimation/fs2000_calibrated_covariance.m.trs
estimation/fs2000_model_comparison.o.trs: estimation/fs2000.o.trs estimation/fs2000_initialize_from_calib.o.trs estimation/fs2000_calibrated_covariance.o.trs
estimation/fs2000_init_from_previous.m.trs: estimation/fs2000.m.trs
estimation/fs2000_init_from_previous.o.trs: estimation/fs2000.o.trs
optimizers/fs2000_102.m.trs: estimation/fs2000.m.trs
optimizers/fs2000_102.o.trs: estimation/fs2000.o.trs

92
tests/estimation/fs2000_init_from_previous.mod Normal file
View File

@ -0,0 +1,92 @@
// Test the mh_initialize_from_previous_mcmc option of the estimation command
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m e_b;
parameters alp bet gam mst rho psi del;
alp = 0.33;
bet = 0.99;
gam = 0.003;
mst = 1.011;
rho = 0.7;
psi = 0.787;
del = 0.02;
model;
dA = exp(gam+e_a);
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
-P/(c(+1)*P(+1)*m)+bet*exp(e_b)*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
W = l/n;
-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
1/(c*P)-bet*exp(e_b)*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
P*c = m;
m-1+d = l;
e = exp(e_a);
y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
gy_obs = dA*y/y(-1);
gp_obs = (P/P(-1))*m(-1)/dA;
end;
steady_state_model;
dA = exp(gam);
gst = 1/dA;
m = mst;
khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
n = xist/(nust+xist);
P = xist + nust;
k = khst*n;
l = psi*mst*n/( (1-psi)*(1-n) );
c = mst/P;
d = l - mst + 1;
y = k^alp*n^(1-alp)*gst^alp;
R = mst/bet;
W = l/n;
ist = y-c;
q = 1 - d;
e = 1;
gp_obs = m/dA;
gy_obs = dA;
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
steady;
check;
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.223;
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_b, inv_gamma_pdf, 0.01, inf;
stderr e_m, inv_gamma_pdf, 0.008862, inf;
end;
varobs gp_obs gy_obs;
% Test mh_initialize_from_previous_mcmc option
estimation(order=1, datafile=fsdat_simul, nobs=192, loglinear, mh_replic=50, mh_nblocks=1, posterior_sampling_method='slice',
mh_initialize_from_previous_mcmc, mh_initialize_from_previous_mcmc_directory='./fs2000');
% Provide record file and prior file
estimation(order=1, datafile=fsdat_simul, nobs=192, loglinear, mh_replic=50, mh_nblocks=1, posterior_sampling_method='slice',
mh_initialize_from_previous_mcmc, mh_initialize_from_previous_mcmc_directory='',
mh_initialize_from_previous_mcmc_record = './fs2000/metropolis/fs2000_mh_history_0.mat',
mh_initialize_from_previous_mcmc_prior = './fs2000/prior/definition.mat');