[WIP] Add Sequential Monte Carlo samplers.
parent
6938080951
commit
a4c77a41e9
GPG Key ID:
295C1FE89E17EB3C
|
|
@ -7271,6 +7271,18 @@ observed variables.
|
|||
Chain draws than the MH-algorithm. Its relative (in)efficiency can be investigated via
|
||||
the reported inefficiency factors.
|
||||
|
||||
``'hssmc'``
|
||||
|
||||
Instructs Dynare to use the *Herbst and Schorfheide (2014)*
|
||||
version of the Sequential Monte-Carlo sampler instead of the
|
||||
standard Random-Walk Metropolis-Hastings.
|
||||
|
||||
``'dsmh'``
|
||||
|
||||
Instructs Dynare to use the Dynamic Striated Metropolis Hastings
|
||||
sampler proposed by *Waggoner, Wu and Zha (2016)* instead of the
|
||||
standard Random-Walk Metropolis-Hastings.
|
||||
|
||||
.. option:: posterior_sampler_options = (NAME, VALUE, ...)
|
||||
|
||||
A list of NAME and VALUE pairs. Can be used to set options for
|
||||
|
|
|
|||
|
|
@ -1,23 +1,22 @@
|
|||
function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, options_, bounds, bayestopt_)
|
||||
|
||||
% function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, options_, bounds, bayestopt_)
|
||||
% initialization of posterior samplers
|
||||
% Initialization of posterior samplers
|
||||
%
|
||||
% INPUTS
|
||||
% posterior_sampler_options: posterior sampler options
|
||||
% options_: structure storing the options
|
||||
% bounds: structure containing prior bounds
|
||||
% bayestopt_: structure storing information about priors
|
||||
|
||||
% - posterior_sampler_options [struct] posterior sampler options
|
||||
% - options_ [struct] options
|
||||
% - bounds [struct] prior bounds
|
||||
% - bayestopt_ [struct] information about priors
|
||||
%
|
||||
% OUTPUTS
|
||||
% posterior_sampler_options: checked posterior sampler options
|
||||
% options_: structure storing the options
|
||||
% bayestopt_: structure storing information about priors
|
||||
% - posterior_sampler_options [struct] checked posterior sampler options (updated)
|
||||
% - options_ [struct] options (updated)
|
||||
% - bayestopt_ [struct] information about priors (updated)
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% none
|
||||
|
||||
% Copyright © 2015-2022 Dynare Team
|
||||
% Copyright © 2015-2023 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
|
|
@ -35,9 +34,9 @@ function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sam
|
|||
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
||||
|
||||
|
||||
init=0;
|
||||
init = false;
|
||||
if isempty(posterior_sampler_options)
|
||||
init=1;
|
||||
init = true;
|
||||
end
|
||||
|
||||
if init
|
||||
|
|
@ -222,7 +221,6 @@ if init
|
|||
end
|
||||
end
|
||||
|
||||
|
||||
case 'slice'
|
||||
posterior_sampler_options.parallel_bar_refresh_rate=1;
|
||||
posterior_sampler_options.serial_bar_refresh_rate=1;
|
||||
|
|
@ -371,6 +369,90 @@ if init
|
|||
posterior_sampler_options.WR=[];
|
||||
end
|
||||
|
||||
case 'hssmc'
|
||||
|
||||
posterior_sampler_options.parallel_bar_refresh_rate=1;
|
||||
posterior_sampler_options.serial_bar_title='HS Sequential Monte-Carlo';
|
||||
|
||||
% default options
|
||||
posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.hssmc);
|
||||
|
||||
% user defined options
|
||||
if ~isempty(options_.posterior_sampler_options.sampling_opt)
|
||||
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
|
||||
for i=1:rows(options_list)
|
||||
switch options_list{i,1}
|
||||
case 'proposal_distribution'
|
||||
if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
|
||||
strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
|
||||
error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
|
||||
'rand_multivariate_student or rand_multivariate_normal as options']);
|
||||
else
|
||||
posterior_sampler_options.proposal_distribution=options_list{i,2};
|
||||
end
|
||||
case 'student_degrees_of_freedom'
|
||||
if options_list{i,2} <= 0
|
||||
error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
|
||||
else
|
||||
posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
|
||||
end
|
||||
case 'save_tmp_file'
|
||||
posterior_sampler_options.save_tmp_file = options_list{i,2};
|
||||
case 'number_of_particles'
|
||||
posterior_sampler_options.nparticles = options_list{i,2};
|
||||
otherwise
|
||||
warning(['rwmh_sampler: Unknown option (' options_list{i,1} ')!'])
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
options_.mode_compute = 0;
|
||||
options_.cova_compute = 0;
|
||||
options_.mh_replic = 0;
|
||||
options_.mh_posterior_mode_estimation = true;
|
||||
|
||||
case 'dsmh'
|
||||
|
||||
posterior_sampler_options.parallel_bar_refresh_rate=1;
|
||||
posterior_sampler_options.serial_bar_title='Dynamic Striated Metropolis Hastings';
|
||||
|
||||
% default options
|
||||
posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.dsmh);
|
||||
|
||||
% user defined options
|
||||
if ~isempty(options_.posterior_sampler_options.sampling_opt)
|
||||
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
|
||||
for i=1:rows(options_list)
|
||||
switch options_list{i,1}
|
||||
case 'proposal_distribution'
|
||||
if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
|
||||
strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
|
||||
error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
|
||||
'rand_multivariate_student or rand_multivariate_normal as options']);
|
||||
else
|
||||
posterior_sampler_options.proposal_distribution=options_list{i,2};
|
||||
end
|
||||
case 'student_degrees_of_freedom'
|
||||
if options_list{i,2} <= 0
|
||||
error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
|
||||
else
|
||||
posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
|
||||
end
|
||||
case 'save_tmp_file'
|
||||
posterior_sampler_options.save_tmp_file = options_list{i,2};
|
||||
case 'number_of_particles'
|
||||
posterior_sampler_options.nparticles = options_list{i,2};
|
||||
otherwise
|
||||
warning(['rwmh_sampler: Unknown option (' options_list{i,1} ')!'])
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
options_.mode_compute = 0;
|
||||
options_.cova_compute = 0;
|
||||
options_.mh_replic = 0;
|
||||
options_.mh_posterior_mode_estimation = true;
|
||||
|
||||
otherwise
|
||||
error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method);
|
||||
end
|
||||
|
|
@ -380,7 +462,7 @@ end
|
|||
|
||||
% here are all samplers requiring a proposal distribution
|
||||
if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
|
||||
if ~options_.cova_compute && ~(options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix)
|
||||
if ~options_.cova_compute && ~(options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix)
|
||||
if strcmp('hessian',options_.MCMC_jumping_covariance)
|
||||
skipline()
|
||||
disp('check_posterior_sampler_options:: I cannot start the MCMC because the Hessian of the posterior kernel at the mode was not computed')
|
||||
|
|
|
|||
|
|
@ -477,17 +477,15 @@ options_.posterior_sampler_options.imh.proposal_distribution = 'rand_multivariat
|
|||
options_.posterior_sampler_options.imh.use_mh_covariance_matrix=0;
|
||||
options_.posterior_sampler_options.imh.save_tmp_file=0;
|
||||
% Herbst and Schorfeide SMC Sampler
|
||||
%options_.posterior_sampler = 'Herbst_Schorfheide' ;
|
||||
options_.posterior_sampler_options.HSsmc.nphi= 25 ;
|
||||
options_.posterior_sampler_options.HSsmc.lambda = 2 ;
|
||||
options_.posterior_sampler_options.HSsmc.nparticles = 20000 ;
|
||||
options_.posterior_sampler_options.HSsmc.c = 0.5 ;
|
||||
options_.posterior_sampler_options.HSsmc.acpt = 0.25 ;
|
||||
options_.posterior_sampler_options.HSsmc.trgt = 0.25 ;
|
||||
options_.posterior_sampler_options.HSsmc.option_mutation = 1 ;
|
||||
options_.posterior_sampler_options.HSsmc.alp = .9 ;
|
||||
options_.posterior_sampler_options.hssmc.nphi= 25 ;
|
||||
options_.posterior_sampler_options.hssmc.lambda = 2 ;
|
||||
options_.posterior_sampler_options.hssmc.nparticles = 20000 ;
|
||||
options_.posterior_sampler_options.hssmc.c = 0.5 ;
|
||||
options_.posterior_sampler_options.hssmc.acpt = 0.25 ;
|
||||
options_.posterior_sampler_options.hssmc.trgt = 0.25 ;
|
||||
options_.posterior_sampler_options.hssmc.option_mutation = true ;
|
||||
options_.posterior_sampler_options.hssmc.alp = .9 ;
|
||||
% DSMH: Dynamic Striated Metropolis-Hastings algorithm
|
||||
%options_.posterior_sampler = 'DSMH' ;
|
||||
options_.posterior_sampler_options.dsmh.H = 25 ;
|
||||
options_.posterior_sampler_options.dsmh.N = 20 ;
|
||||
options_.posterior_sampler_options.dsmh.G = 10 ;
|
||||
|
|
@ -495,7 +493,7 @@ options_.posterior_sampler_options.dsmh.K = 50 ;
|
|||
options_.posterior_sampler_options.dsmh.lambda1 = 0.1 ;
|
||||
options_.posterior_sampler_options.dsmh.nparticles = 20000 ;
|
||||
options_.posterior_sampler_options.dsmh.alpha0 = 0.2 ;
|
||||
options_.posterior_sampler_options.dsmh.alpha1 = 0.3 ;
|
||||
options_.posterior_sampler_options.dsmh.alpha1 = 0.3 ;
|
||||
options_.posterior_sampler_options.dsmh.tau = 10 ;
|
||||
|
||||
options_.trace_plot_ma = 200;
|
||||
|
|
|
|||
|
|
@ -1,8 +1,6 @@
|
|||
function [dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_, bounds] = dynare_estimation_init(var_list_, dname, gsa_flag, M_, options_, oo_, estim_params_, bayestopt_)
|
||||
|
||||
% function [dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_, bounds] = dynare_estimation_init(var_list_, dname, gsa_flag, M_, options_, oo_, estim_params_, bayestopt_)
|
||||
% performs initialization tasks before estimation or
|
||||
% global sensitivity analysis
|
||||
% Performs initialization tasks before estimation or global sensitivity analysis
|
||||
%
|
||||
% INPUTS
|
||||
% var_list_: selected endogenous variables vector
|
||||
|
|
@ -531,13 +529,13 @@ end
|
|||
%set options for old interface from the ones for new interface
|
||||
if ~isempty(dataset_)
|
||||
options_.nobs = dataset_.nobs;
|
||||
if options_.endogenous_prior
|
||||
if options_.endogenous_prior
|
||||
if dataset_info.missing.no_more_missing_observations<dataset_.nobs-10
|
||||
fprintf('\ndynare_estimation_init: There are missing observations in the data.\n')
|
||||
fprintf('dynare_estimation_init: I am computing the moments for the endogenous prior only\n')
|
||||
fprintf('dynare_estimation_init: on the observations after the last missing one, i.e. %u.\n',dataset_info.missing.no_more_missing_observations)
|
||||
else
|
||||
fprintf('\ndynare_estimation_init: There are too many missing observations in the data.\n')
|
||||
fprintf('\ndynare_estimation_init: There are too many missing observations in the data.\n')
|
||||
fprintf('dynare_estimation_init: The endogenous_prior-option needs a consistent sample of \n')
|
||||
fprintf('dynare_estimation_init: at least 10 full observations at the end.\n')
|
||||
error('The endogenous_prior-option does not support your missing data.')
|
||||
|
|
@ -676,15 +674,15 @@ end
|
|||
|
||||
if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter
|
||||
if ~isempty(options_.nk)
|
||||
fprintf('dynare_estimation_init: the inversion filter does not support filter_step_ahead. Disabling the option.\n')
|
||||
fprintf('dynare_estimation_init: the inversion filter does not support filter_step_ahead. Disabling the option.\n')
|
||||
options_.nk=[];
|
||||
end
|
||||
if options_.filter_covariance
|
||||
fprintf('dynare_estimation_init: the inversion filter does not support filter_covariance. Disabling the option.\n')
|
||||
fprintf('dynare_estimation_init: the inversion filter does not support filter_covariance. Disabling the option.\n')
|
||||
options_.filter_covariance=false;
|
||||
end
|
||||
if options_.smoothed_state_uncertainty
|
||||
fprintf('dynare_estimation_init: the inversion filter does not support smoothed_state_uncertainty. Disabling the option.\n')
|
||||
fprintf('dynare_estimation_init: the inversion filter does not support smoothed_state_uncertainty. Disabling the option.\n')
|
||||
options_.smoothed_state_uncertainty=false;
|
||||
end
|
||||
end
|
||||
|
|
|
|||
|
|
@ -1,39 +1,23 @@
|
|||
function Herbst_Schorfheide_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
|
||||
% function Herbst_Schorfheide_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
|
||||
% SMC sampler from JAE 2014 .
|
||||
function hssmc(TargetFun, xparam1, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
|
||||
|
||||
% Sequential Monte-Carlo sampler, Herbst and Schorfheide (JAE, 2014).
|
||||
%
|
||||
% INPUTS
|
||||
% o TargetFun [char] string specifying the name of the objective
|
||||
% function (posterior kernel).
|
||||
% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values).
|
||||
% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters.
|
||||
% o dataset_ data structure
|
||||
% o dataset_info dataset info structure
|
||||
% o options_ options structure
|
||||
% o M_ model structure
|
||||
% o estim_params_ estimated parameters structure
|
||||
% o bayestopt_ estimation options structure
|
||||
% o oo_ outputs structure
|
||||
% - TargetFun [char] string specifying the name of the objective function (posterior kernel).
|
||||
% - xparam1 [double] p×1 vector of parameters to be estimated (initial values).
|
||||
% - mh_bounds [double] p×2 matrix defining lower and upper bounds for the parameters.
|
||||
% - dataset_ [dseries] sample
|
||||
% - dataset_info [struct] informations about the dataset
|
||||
% - options_ [struct] dynare's options
|
||||
% - M_ [struct] model description
|
||||
% - estim_params_ [struct] estimated parameters
|
||||
% - bayestopt_ [struct] estimated parameters
|
||||
% - oo_ [struct] outputs
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% None.
|
||||
%
|
||||
% PARALLEL CONTEXT
|
||||
% The most computationally intensive part of this function may be executed
|
||||
% in parallel. The code suitable to be executed in
|
||||
% parallel on multi core or cluster machine (in general a 'for' cycle)
|
||||
% has been removed from this function and been placed in the posterior_sampler_core.m funtion.
|
||||
%
|
||||
% The DYNARE parallel packages comprise a i) set of pairs of Matlab functions that can be executed in
|
||||
% parallel and called name_function.m and name_function_core.m and ii) a second set of functions used
|
||||
% to manage the parallel computations.
|
||||
%
|
||||
% This function was the first function to be parallelized. Later, other
|
||||
% functions have been parallelized using the same methodology.
|
||||
% Then the comments write here can be used for all the other pairs of
|
||||
% parallel functions and also for management functions.
|
||||
% None.
|
||||
|
||||
% Copyright © 2006-2022 Dynare Team
|
||||
% Copyright © 2022-2023 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
|
|
@ -50,17 +34,21 @@ function Herbst_Schorfheide_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset
|
|||
% You should have received a copy of the GNU General Public License
|
||||
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
||||
|
||||
opts = options_.posterior_sampler_options.hssmc;
|
||||
|
||||
% Create the tempering schedule
|
||||
phi = bsxfun(@power,(bsxfun(@minus,1:1:options_.posterior_sampler_options.HSsmc.nphi,1)/(options_.posterior_sampler_options.HSsmc.nphi-1)),options_.posterior_sampler_options.HSsmc.lambda) ;
|
||||
phi = bsxfun(@power, (0:opts.nphi-1)/(opts.nphi-1), opts.lambda);
|
||||
|
||||
% tuning for MH algorithms matrices
|
||||
zhat = 0 ; % normalization constant
|
||||
csim = zeros(options_.posterior_sampler_options.HSsmc.nphi,1) ; % scale parameter
|
||||
ESSsim = zeros(options_.posterior_sampler_options.HSsmc.nphi,1) ; % ESS
|
||||
acptsim = zeros(options_.posterior_sampler_options.HSsmc.nphi,1) ; % average acceptance rate
|
||||
zhat = 0; % normalization constant
|
||||
csim = zeros(opts.nphi,1); % scale parameter
|
||||
ESSsim = zeros(opts.nphi,1); % ESS
|
||||
acptsim = zeros(opts.nphi,1); % average acceptance rate
|
||||
|
||||
% Step 0: Initialization of the sampler
|
||||
[ param, tlogpost_i, loglik, bayestopt_] = ...
|
||||
SMC_samplers_initialization(TargetFun, xparam1, mh_bounds, dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,options_.posterior_sampler_options.HSsmc.nparticles);
|
||||
|
||||
weights = ones(options_.posterior_sampler_options.HSsmc.nparticles,1)/options_.posterior_sampler_options.HSsmc.nparticles ;
|
||||
% The Herbst and Schorfheide sampler starts here
|
||||
for i=2:options_.posterior_sampler_options.HSsmc.nphi
|
||||
|
|
@ -105,6 +93,7 @@ for i=2:options_.posterior_sampler_options.HSsmc.nphi
|
|||
fprintf(' Neff = %5.4f \n', ESSsim(i));
|
||||
fprintf(' %accept. = %5.4f \n', acptsim(i));
|
||||
end
|
||||
|
||||
indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.HSsmc.nparticles);
|
||||
distrib_param = param(:,indx_resmpl);
|
||||
fprintf(' Log_lik = %5.4f \n', zhat);
|
||||
|
|
@ -159,7 +148,7 @@ for k=1:npar %min(nstar,npar-(plt-1)*nstar)
|
|||
[name,texname] = get_the_name(k,TeX,M_,estim_params_,options_);
|
||||
optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.posterior_sampler_options.HSsmc.nparticles,bandwidth,kernel_function);
|
||||
[density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
|
||||
options_.posterior_sampler_options.HSsmc.nparticles,optimal_bandwidth,kernel_function);
|
||||
options_.posterior_sampler_options.HSsmc.nparticles,optimal_bandwidth,kernel_function);
|
||||
plot(density(:,1),density(:,2));
|
||||
hold on
|
||||
if TeX
|
||||
|
|
@ -1,8 +1,7 @@
|
|||
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, bayestopt_] = ...
|
||||
% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
|
||||
% Metropolis-Hastings initialization.
|
||||
|
||||
% Posterior sampler initialization.
|
||||
%
|
||||
% INPUTS
|
||||
% o TargetFun [char] string specifying the name of the objective
|
||||
|
|
@ -132,7 +131,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
|
|||
MetropolisFolder0 = fileparts(RecordFile0);
|
||||
PreviousFolder0=fileparts(MetropolisFolder0);
|
||||
[~, ModelName0]=fileparts(PreviousFolder0);
|
||||
else
|
||||
else
|
||||
%% check for proper filesep char in user defined paths
|
||||
PreviousFolder0=strrep(PreviousFolder0,'\',filesep);
|
||||
MetropolisFolder0 = [PreviousFolder0 filesep 'metropolis'];
|
||||
|
|
@ -228,7 +227,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
|
|||
if options_.nointeractive
|
||||
disp(['Estimation::mcmc: I reduce mh_init_scale_factor by 10 percent:'])
|
||||
if isfield(options_,'mh_init_scale')
|
||||
options_.mh_init_scale = .9*options_.mh_init_scale;
|
||||
options_.mh_init_scale = .9*options_.mh_init_scale;
|
||||
fprintf('Estimation::mcmc: Parameter mh_init_scale is now equal to %f.\n',options_.mh_init_scale)
|
||||
else
|
||||
options_.mh_init_scale_factor = .9*options_.mh_init_scale_factor;
|
||||
|
|
@ -301,7 +300,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
|
|||
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);
|
||||
|
|
@ -502,7 +501,7 @@ elseif options_.mh_recover
|
|||
% How many mh-files are saved in this block?
|
||||
ExistingDrawsInLastMCFile=zeros(NumberOfBlocks,1); %initialize: no MCMC draws of current MCMC are in file from last run
|
||||
% Check whether last present file is a file included in the last MCMC run
|
||||
|
||||
|
||||
update_record=0;
|
||||
for k=1:NumberOfBlocks
|
||||
FirstBlock = k;
|
||||
|
|
|
|||
|
|
@ -0,0 +1,85 @@
|
|||
// See fs2000.mod in the examples/ directory for details on the model
|
||||
|
||||
var m P c e W R k d n l gy_obs gp_obs y dA;
|
||||
varexo e_a e_m;
|
||||
|
||||
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*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*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_m, inv_gamma_pdf, 0.008862, inf;
|
||||
end;
|
||||
|
||||
varobs gp_obs gy_obs;
|
||||
|
||||
options_.solve_tolf = 1e-12;
|
||||
|
||||
estimation(order=1,datafile='../fsdat_simul.m',nobs=192,loglinear,posterior_sampling_method='dsmh');
|
||||
|
|
@ -0,0 +1,85 @@
|
|||
// See fs2000.mod in the examples/ directory for details on the model
|
||||
|
||||
var m P c e W R k d n l gy_obs gp_obs y dA;
|
||||
varexo e_a e_m;
|
||||
|
||||
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*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*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_m, inv_gamma_pdf, 0.008862, inf;
|
||||
end;
|
||||
|
||||
varobs gp_obs gy_obs;
|
||||
|
||||
options_.solve_tolf = 1e-12;
|
||||
|
||||
estimation(order=1,datafile='../fsdat_simul.m',nobs=192,loglinear,posterior_sampling_method='hssmc');
|
||||
Loading…
Reference in New Issue