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

Add new class for priors. 

This commit only changes the routine used to draw random deviates from
the prior distribution, without relying on persistent variables (which allows
parallelisation).
mr#2134
Stéphane Adjemian (Ryûk) 2023-04-26 10:34:25 +02:00
parent 339ba7102e
commit 015513380f
Signed by: stepan
GPG Key ID: 295C1FE89E17EB3C
12 changed files with 485 additions and 82 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

376
matlab/@dprior/dprior.m Normal file
View File

@ -0,0 +1,376 @@
classdef dprior
properties
p6 = []; % Prior first hyperparameter.
p7 = []; % Prior second hyperparameter.
p3 = []; % Lower bound of the prior support.
p4 = []; % Upper bound of the prior support.
lb = []; % Truncated prior lower bound.
ub = []; % Truncated prior upper bound.
uniform_index = []; % Index for the uniform priors.
gaussian_index = []; % Index for the gaussian priors.
gamma_index = []; % Index for the gamma priors.
beta_index = []; % Index for the beta priors.
inverse_gamma_1_index = []; % Index for the inverse gamma type 1 priors.
inverse_gamma_2_index = []; % Index for the inverse gamma type 2 priors.
weibull_index = []; % Index for the weibull priors.
uniform_draws = false;
gaussian_draws = false;
gamma_draws = false;
beta_draws = false;
inverse_gamma_1_draws = false;
inverse_gamma_2_draws = false;
weibull_draws = false;
end
methods
function o = dprior(BayesInfo, PriorTrunc, Uniform)
% Class constructor.
%
% INPUTS
% - BayesInfo [struct] Informations about the prior distribution, aka bayestopt_.
% - PriorTrunc [double] scalar, probability mass to be excluded, aka options_.prior_trunc
% - Uniform [logical] scalar, produce uniform random deviates on the prior support.
%
% OUTPUTS
% - o [dprior] scalar, prior object.
%
% REQUIREMENTS
% None.
o.p6 = BayesInfo.p6;
o.p7 = BayesInfo.p7;
o.p3 = BayesInfo.p3;
o.p4 = BayesInfo.p4;
bounds = prior_bounds(BayesInfo, PriorTrunc);
o.lb = bounds.lb;
o.ub = bounds.ub;
if nargin>2 && Uniform
prior_shape = repmat(5, length(o.p6), 1);
else
prior_shape = BayesInfo.pshape;
end
o.beta_index = find(prior_shape==1);
if ~isempty(o.beta_index)
o.beta_draws = true;
end
o.gamma_index = find(prior_shape==2);
if ~isempty(o.gamma_index)
o.gamma_draws = true;
end
o.gaussian_index = find(prior_shape==3);
if ~isempty(o.gaussian_index)
o.gaussian_draws = true;
end
o.inverse_gamma_1_index = find(prior_shape==4);
if ~isempty(o.inverse_gamma_1_index)
o.inverse_gamma_1_draws = true;
end
o.uniform_index = find(prior_shape==5);
if ~isempty(o.uniform_index)
o.uniform_draws = true;
end
o.inverse_gamma_2_index = find(prior_shape==6);
if ~isempty(o.inverse_gamma_2_index)
o.inverse_gamma_2_draws = true;
end
o.weibull_index = find(prior_shape==8);
if ~isempty(o.weibull_index)
o.weibull_draws = true;
end
end
function p = draw(o)
% Return a random draw from the prior distribution.
%
% INPUTS
% - o [dprior]
%
% OUTPUTS
% - p [double] m×1 vector, random draw from the prior distribution (m is the number of estimated parameters).
%
% REMARKS
% None.
%
% EXAMPLE
%
% >> Prior = dprior(bayestopt_, options_.prior_trunc);
% >> d = Prior.draw()
p = NaN(rows(o.lb), 1);
if o.uniform_draws
p(o.uniform_index) = rand(length(o.uniform_index),1).*(o.p4(o.uniform_index)-o.p3(o.uniform_index)) + o.p3(o.uniform_index);
out_of_bound = find( (p(o.uniform_index)>o.ub(o.uniform_index)) | (p(o.uniform_index)<o.lb(o.uniform_index)));
while ~isempty(out_of_bound)
p(o.uniform_index) = rand(length(o.uniform_index), 1).*(o.p4(o.uniform_index)-o.p3(o.uniform_index)) + o.p3(o.uniform_index);
out_of_bound = find( (p(o.uniform_index)>o.ub(o.uniform_index)) | (p(o.uniform_index)<o.lb(o.uniform_index)));
end
end
if o.gaussian_draws
p(o.gaussian_index) = randn(length(o.gaussian_index), 1).*o.p7(o.gaussian_index) + o.p6(o.gaussian_index);
out_of_bound = find( (p(o.gaussian_index)>o.ub(o.gaussian_index)) | (p(o.gaussian_index)<o.lb(o.gaussian_index)));
while ~isempty(out_of_bound)
p(o.gaussian_index(out_of_bound)) = randn(length(o.gaussian_index(out_of_bound)), 1).*o.p7(o.gaussian_index(out_of_bound)) + o.p6(o.gaussian_index(out_of_bound));
out_of_bound = find( (p(o.gaussian_index)>o.ub(o.gaussian_index)) | (p(o.gaussian_index)<o.lb(o.gaussian_index)));
end
end
if o.gamma_draws
p(o.gamma_index) = gamrnd(o.p6(o.gamma_index), o.p7(o.gamma_index))+o.p3(o.gamma_index);
out_of_bound = find( (p(o.gamma_index)>o.ub(o.gamma_index)) | (p(o.gamma_index)<o.lb(o.gamma_index)));
while ~isempty(out_of_bound)
p(o.gamma_index(out_of_bound)) = gamrnd(o.p6(o.gamma_index(out_of_bound)), o.p7(o.gamma_index(out_of_bound)))+o.p3(o.gamma_index(out_of_bound));
out_of_bound = find( (p(o.gamma_index)>o.ub(o.gamma_index)) | (p(o.gamma_index)<o.lb(o.gamma_index)));
end
end
if o.beta_draws
p(o.beta_index) = (o.p4(o.beta_index)-o.p3(o.beta_index)).*betarnd(o.p6(o.beta_index), o.p7(o.beta_index))+o.p3(o.beta_index);
out_of_bound = find( (p(o.beta_index)>o.ub(o.beta_index)) | (p(o.beta_index)<o.lb(o.beta_index)));
while ~isempty(out_of_bound)
p(o.beta_index(out_of_bound)) = (o.p4(o.beta_index(out_of_bound))-o.p3(o.beta_index(out_of_bound))).*betarnd(o.p6(o.beta_index(out_of_bound)), o.p7(o.beta_index(out_of_bound)))+o.p3(o.beta_index(out_of_bound));
out_of_bound = find( (p(o.beta_index)>o.ub(o.beta_index)) | (p(o.beta_index)<o.lb(o.beta_index)));
end
end
if o.inverse_gamma_1_draws
p(o.inverse_gamma_1_index) = ...
sqrt(1./gamrnd(o.p7(o.inverse_gamma_1_index)/2, 2./o.p6(o.inverse_gamma_1_index)))+o.p3(o.inverse_gamma_1_index);
out_of_bound = find( (p(o.inverse_gamma_1_index)>o.ub(o.inverse_gamma_1_index)) | (p(o.inverse_gamma_1_index)<o.lb(o.inverse_gamma_1_index)));
while ~isempty(out_of_bound)
p(o.inverse_gamma_1_index(out_of_bound)) = ...
sqrt(1./gamrnd(o.p7(o.inverse_gamma_1_index(out_of_bound))/2, 2./o.p6(o.inverse_gamma_1_index(out_of_bound))))+o.p3(o.inverse_gamma_1_index(out_of_bound));
out_of_bound = find( (p(o.inverse_gamma_1_index)>o.ub(inverse_gamma_1_index)) | (p(o.inverse_gamma_1_index)<o.lb(o.inverse_gamma_1_index)));
end
end
if o.inverse_gamma_2_draws
p(o.inverse_gamma_2_index) = ...
1./gamrnd(o.p7(o.inverse_gamma_2_index)/2, 2./o.p6(o.inverse_gamma_2_index))+o.p3(o.inverse_gamma_2_index);
out_of_bound = find( (p(o.inverse_gamma_2_index)>o.ub(o.inverse_gamma_2_index)) | (p(o.inverse_gamma_2_index)<o.lb(o.inverse_gamma_2_index)));
while ~isempty(out_of_bound)
p(o.inverse_gamma_2_index(out_of_bound)) = ...
1./gamrnd(o.p7(o.inverse_gamma_2_index(out_of_bound))/2, 2./o.p6(o.inverse_gamma_2_index(out_of_bound)))+o.p3(o.inverse_gamma_2_index(out_of_bound));
out_of_bound = find( (p(o.inverse_gamma_2_index)>o.ub(o.inverse_gamma_2_index)) | (p(o.inverse_gamma_2_index)<o.lb(o.inverse_gamma_2_index)));
end
end
if o.weibull_draws
p(o.weibull_index) = wblrnd(o.p7(o.weibull_index), o.p6(o.weibull_index)) + o.p3(o.weibull_index);
out_of_bound = find( (p(o.weibull_index)>o.ub(o.weibull_index)) | (p(o.weibull_index)<o.lb(o.weibull_index)));
while ~isempty(out_of_bound)
p(o.weibull_index(out_of_bound)) = wblrnd(o.p7(o.weibull_index(out_of_bound)), o.p6(o.weibull_index(out_of_bound)))+o.p3(o.weibull_index(out_of_bound));
out_of_bound = find( (p(o.weibull_index)>o.ub(o.weibull_index)) | (p(o.weibull_index)<o.lb(o.weibull_index)));
end
end
end
function P = draws(o, n)
% Return n independent random draws from the prior distribution.
%
% INPUTS
% - o [dprior]
%
% OUTPUTS
% - P [double] m×n matrix, random draw from the prior distribution.
%
% REMARKS
% If the Parallel Computing Toolbox is available, the main loop is run in parallel.
%
% EXAMPLE
%
% >> Prior = dprior(bayestopt_, options_.prior_trunc);
% >> Prior.draws(1e6)
P = NaN(rows(o.lb), 1);
parfor i=1:n
P(:,i) = draw(o);
end
end
end % methods
end % classdef --*-- Unit tests --*--
%@test:1
%$ % Fill global structures with required fields...
%$ prior_trunc = 1e-10;
%$ p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
%$ p1 = .4*ones(14,1); % Prior mean
%$ p2 = .2*ones(14,1); % Prior std.
%$ p3 = NaN(14,1);
%$ p4 = NaN(14,1);
%$ p5 = NaN(14,1);
%$ p6 = NaN(14,1);
%$ p7 = NaN(14,1);
%$
%$ for i=1:14
%$ switch p0(i)
%$ case 1
%$ % Beta distribution
%$ p3(i) = 0;
%$ p4(i) = 1;
%$ [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
%$ case 2
%$ % Gamma distribution
%$ p3(i) = 0;
%$ p4(i) = Inf;
%$ [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
%$ case 3
%$ % Normal distribution
%$ p3(i) = -Inf;
%$ p4(i) = Inf;
%$ p6(i) = p1(i);
%$ p7(i) = p2(i);
%$ p5(i) = p1(i);
%$ case 4
%$ % Inverse Gamma (type I) distribution
%$ p3(i) = 0;
%$ p4(i) = Inf;
%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
%$ case 5
%$ % Uniform distribution
%$ [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
%$ p3(i) = p6(i);
%$ p4(i) = p7(i);
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
%$ case 6
%$ % Inverse Gamma (type II) distribution
%$ p3(i) = 0;
%$ p4(i) = Inf;
%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
%$ case 8
%$ % Weibull distribution
%$ p3(i) = 0;
%$ p4(i) = Inf;
%$ [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
%$ otherwise
%$ error('This density is not implemented!')
%$ end
%$ end
%$
%$ BayesInfo.pshape = p0;
%$ BayesInfo.p1 = p1;
%$ BayesInfo.p2 = p2;
%$ BayesInfo.p3 = p3;
%$ BayesInfo.p4 = p4;
%$ BayesInfo.p5 = p5;
%$ BayesInfo.p6 = p6;
%$ BayesInfo.p7 = p7;
%$
%$ ndraws = 1e5;
%$ m0 = BayesInfo.p1; %zeros(14,1);
%$ v0 = diag(BayesInfo.p2.^2); %zeros(14);
%$
%$ % Call the tested routine
%$ try
%$ % Instantiate dprior object
%$ o = dprior(BayesInfo, prior_trunc, false);
%$ % Do simulations in a loop and estimate recursively the mean and the variance.
%$ for i = 1:ndraws
%$ draw = o.draw();
%$ m1 = m0 + (draw-m0)/i;
%$ m2 = m1*m1';
%$ v0 = v0 + ((draw*draw'-m2-v0) + (i-1)*(m0*m0'-m2'))/i;
%$ m0 = m1;
%$ end
%$ t(1) = true;
%$ catch
%$ t(1) = false;
%$ end
%$
%$ if t(1)
%$ t(2) = all(abs(m0-BayesInfo.p1)<3e-3);
%$ t(3) = all(all(abs(v0-diag(BayesInfo.p2.^2))<5e-3));
%$ end
%$ T = all(t);
%@eof:1
%@test:2
%$ % Fill global structures with required fields...
%$ prior_trunc = 1e-10;
%$ p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
%$ p1 = .4*ones(14,1); % Prior mean
%$ p2 = .2*ones(14,1); % Prior std.
%$ p3 = NaN(14,1);
%$ p4 = NaN(14,1);
%$ p5 = NaN(14,1);
%$ p6 = NaN(14,1);
%$ p7 = NaN(14,1);
%$
%$ for i=1:14
%$ switch p0(i)
%$ case 1
%$ % Beta distribution
%$ p3(i) = 0;
%$ p4(i) = 1;
%$ [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
%$ case 2
%$ % Gamma distribution
%$ p3(i) = 0;
%$ p4(i) = Inf;
%$ [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
%$ case 3
%$ % Normal distribution
%$ p3(i) = -Inf;
%$ p4(i) = Inf;
%$ p6(i) = p1(i);
%$ p7(i) = p2(i);
%$ p5(i) = p1(i);
%$ case 4
%$ % Inverse Gamma (type I) distribution
%$ p3(i) = 0;
%$ p4(i) = Inf;
%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
%$ case 5
%$ % Uniform distribution
%$ [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
%$ p3(i) = p6(i);
%$ p4(i) = p7(i);
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
%$ case 6
%$ % Inverse Gamma (type II) distribution
%$ p3(i) = 0;
%$ p4(i) = Inf;
%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
%$ case 8
%$ % Weibull distribution
%$ p3(i) = 0;
%$ p4(i) = Inf;
%$ [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
%$ otherwise
%$ error('This density is not implemented!')
%$ end
%$ end
%$
%$ BayesInfo.pshape = p0;
%$ BayesInfo.p1 = p1;
%$ BayesInfo.p2 = p2;
%$ BayesInfo.p3 = p3;
%$ BayesInfo.p4 = p4;
%$ BayesInfo.p5 = p5;
%$ BayesInfo.p6 = p6;
%$ BayesInfo.p7 = p7;
%$
%$ ndraws = 1e5;
%$
%$ % Call the tested routine
%$ try
%$ % Instantiate dprior object.
%$ o = dprior(BayesInfo, prior_trunc, false);
%$ X = o.draws(ndraws);
%$ m = mean(X, 2);
%$ v = var(X, 0, 2);
%$ t(1) = true;
%$ catch
%$ t(1) = false;
%$ end
%$
%$ if t(1)
%$ t(2) = all(abs(m-BayesInfo.p1)<3e-3);
%$ t(3) = all(all(abs(v-BayesInfo.p2.^2)<5e-3));
%$ end
%$ T = all(t);
%@eof:2

21
matlab/GetOneDraw.m
View File

@ -1,5 +1,5 @@
function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_)
% function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_)
function [xparams, logpost] = GetOneDraw(distribution, M_, estim_params_, oo_, options_, bayestopt_)
% draws one parameter vector and its posterior from MCMC or the prior
%
% INPUTS
@ -18,7 +18,7 @@ function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,baye
% SPECIAL REQUIREMENTS
% none
% Copyright © 2005-2017 Dynare Team
% Copyright © 2005-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -35,12 +35,13 @@ function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,baye
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
switch type
case 'posterior'
[xparams, logpost] = metropolis_draw(0);
case 'prior'
xparams = prior_draw();
if isprior(distribution)
xparams = distribution.draw();
if nargout>1
logpost = evaluate_posterior_kernel(xparams',M_,estim_params_,oo_,options_,bayestopt_);
logpost = evaluate_posterior_kernel(xparams, M_, estim_params_, oo_, options_, bayestopt_);
end
end
elseif ischar(distribution) && strcmpi(distribution, 'posterior')
[xparams, logpost] = metropolis_draw(0);
else
error('GetOneDraw:: Wrong inputs.')
end

13
matlab/PosteriorIRF_core1.m
View File

@ -23,7 +23,7 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab)
% SPECIAL REQUIREMENTS.
% None.
%
% Copyright © 2006-2019 Dynare Team
% Copyright © 2006-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -134,12 +134,17 @@ end
% Parallel 'while' very good!!!
stock_param=zeros(MAX_nruns,npar);
stock_irf_dsge=zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge);
if strcmp(type, 'prior')
Prior = dprior(bayestopt_, options_.prior_trunc);
end
while fpar<B
fpar = fpar + 1;
irun = irun+1;
irun2 = irun2+1;
if strcmpi(type,'prior')
deep = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
deep = Prior.draw();
else
deep = x(fpar,:);
end
@ -294,7 +299,7 @@ end
% directory on call machine that contain the model).
myoutput.OutputFileName = [OutputFileName_dsge;
OutputFileName_param;
OutputFileName_bvardsge];
OutputFileName_param;
OutputFileName_bvardsge];
myoutput.nosaddle = nosaddle;

21
matlab/dynare_identification.m
View File

@ -40,13 +40,14 @@ function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, ST
% * identification_analysis
% * isoctave
% * plot_identification
% * prior_draw
% * dprior.draw
% * set_default_option
% * set_prior
% * skipline
% * vnorm
% =========================================================================
% Copyright © 2010-2022 Dynare Team
% Copyright © 2010-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -328,18 +329,18 @@ if options_.discretionary_policy || options_.ramsey_policy
options_ident.analytic_derivation_mode=-1;
end
end
options_.analytic_derivation_mode = options_ident.analytic_derivation_mode; %overwrite setting in options_
% initialize persistent variables in prior_draw
% Instantiate dprior object (Prior)
if prior_exist
if any(bayestopt_.pshape > 0)
if options_ident.prior_range
%sample uniformly from prior ranges (overwrite prior specification)
prior_draw(bayestopt_, options_.prior_trunc, true);
Prior = dprior(bayestopt_, options_.prior_trunc, true);
else
%sample from prior distributions
prior_draw(bayestopt_, options_.prior_trunc, false);
Prior = dprior(bayestopt_, options_.prior_trunc, false);
end
else
options_ident.prior_mc = 1; %only one single point
@ -495,7 +496,7 @@ if iload <=0
kk=0;
while kk<50 && info(1)
kk=kk+1;
params = prior_draw();
params = Prior.draw();
options_ident.tittxt = 'Random_prior_params'; %title text for graphs and figures
% perform identification analysis
[ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, error_indicator_point] = ...
@ -548,7 +549,7 @@ if iload <=0
if external_sample
params = pdraws0(iteration+1,:); % loaded draws
else
params = prior_draw(); % new random draw from prior
params = Prior.draw(); % new random draw from prior
end
options_ident.tittxt = []; % clear title text for graphs and figures
% run identification analysis
@ -692,7 +693,7 @@ if iload <=0
end
% store results for minimal system
if ~options_MC.no_identification_minimal
if ~options_MC.no_identification_minimal
if ~error_indicator.identification_minimal
STO_dMINIMAL(:,:,run_index) = ide_minimal.dMINIMAL;
IDE_MINIMAL.cond(iteration,1) = ide_minimal.cond;
@ -901,7 +902,7 @@ if SampleSize > 1
store_nodisplay = options_.nodisplay;
options_.nodisplay = 1;
% HIGHEST CONDITION NUMBER
[~, jmax] = max(IDE_MOMENTS.cond);
[~, jmax] = max(IDE_MOMENTS.cond);
tittxt = 'Draw with HIGHEST condition number';
fprintf('\nTesting %s.\n',tittxt);
if ~iload

56
matlab/execute_prior_posterior_function.m
View File

@ -17,7 +17,7 @@ function oo_=execute_prior_posterior_function(posterior_function_name,M_,options
% OUTPUTS
% oo_ [structure] Matlab/Octave structure gathering the results (initialized by dynare, see @ref{oo_}).
% Copyright © 2013-2015 Dynare Team
% Copyright © 2013-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -47,21 +47,21 @@ end
%Create function handle
functionhandle=str2func(posterior_function_name);
prior = true;
n_draws=options_.sampling_draws;
% Get informations about the _posterior_draws files.
if strcmpi(type,'posterior')
%% discard first mh_drop percent of the draws:
% Get informations about the _posterior_draws files.
% discard first mh_drop percent of the draws:
CutSample(M_, options_, estim_params_);
%% initialize metropolis draws
options_.sub_draws=n_draws; %set draws for sampling; changed value is not returned to base workspace
[error_flag,~,options_]= metropolis_draw(1,options_,estim_params_,M_);
% initialize metropolis draws
options_.sub_draws = n_draws; % set draws for sampling; changed value is not returned to base workspace
[error_flag, ~, options_] = metropolis_draw(1, options_, estim_params_, M_);
if error_flag
error('EXECUTE_POSTERIOR_FUNCTION: The draws could not be initialized')
end
n_draws=options_.sub_draws;
prior = false;
n_draws = options_.sub_draws;
elseif strcmpi(type,'prior')
% Get informations about the prior distribution.
if isempty(bayestopt_)
if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && (size(estim_params_.var_exo,1)+size(estim_params_.var_endo,1)+size(estim_params_.corrx,1)+size(estim_params_.corrn,1)+size(estim_params_.param_vals,1))==0)
[xparam1,estim_params_,bayestopt_,lb,ub,M_] = set_prior(estim_params_,M_,options_);
@ -72,38 +72,36 @@ elseif strcmpi(type,'prior')
if exist([M_.fname '_prior_restrictions.m'])
warning('prior_function currently does not support endogenous prior restrictions. They will be ignored. Consider using a prior_function with nobs=1.')
end
prior_draw(bayestopt_, options_.prior_trunc);
Prior = dprior(bayestopt_, options_.prior_trunc);
else
error('EXECUTE_POSTERIOR_FUNCTION: Unknown type!')
end
%get draws for later use
first_draw=GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
parameter_mat=NaN(n_draws,length(first_draw));
parameter_mat(1,:)=first_draw;
for draw_iter=2:n_draws
parameter_mat(draw_iter,:) = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
if strcmpi(type, 'prior')
parameter_mat = Prior.draws(n_draws);
else
parameter_mat = NaN(length(bayestopt_.p6), n_draws);
for i = 1:n_draws
parameter_mat(:,i) = GetOneDraw(type, M_, estim_params_, oo_, options_, bayestopt_);
end
end
% get output size
% Get output size
try
junk=functionhandle(parameter_mat(1,:),M_,options_,oo_,estim_params_,bayestopt_,dataset_,dataset_info);
junk = functionhandle(parameter_mat(:,1), M_, options_, oo_, estim_params_, bayestopt_, dataset_, dataset_info);
catch err
fprintf('\nEXECUTE_POSTERIOR_FUNCTION: Execution of prior/posterior function led to an error. Execution cancelled.\n')
rethrow(err)
end
%initialize cell with number of columns
results_cell=cell(n_draws,size(junk,2));
% Initialize cell with number of columns
results_cell = cell(n_draws, columns(junk));
%% compute function on draws
for draw_iter = 1:n_draws
M_ = set_all_parameters(parameter_mat(draw_iter,:),estim_params_,M_);
[results_cell(draw_iter,:)]=functionhandle(parameter_mat(draw_iter,:),M_,options_,oo_,estim_params_,bayestopt_,dataset_,dataset_info);
% Evaluate function on each draw
for i = 1:n_draws
M_ = set_all_parameters(parameter_mat(:,i), estim_params_, M_);
[results_cell(i,:)] = functionhandle(parameter_mat(:,i), M_, options_, oo_, estim_params_, bayestopt_, dataset_, dataset_info);
end
if prior
oo_.prior_function_results = results_cell;
else
oo_.posterior_function_results = results_cell;
end
% Save results under oo_
oo_.(sprintf('%s_function_results', type)) = results_cell;

22
matlab/isprior.m Normal file
View File

@ -0,0 +1,22 @@
function b = isprior(o)
% Return True iff o is a dprior object.
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
b = isa(o, 'dprior');

2
matlab/particles

@ -1 +1 @@
Subproject commit 7926e75aa525b73c91bbb944d19556c43a8fbbe6
Subproject commit bb16df7f2e6ddbc81e23b9661ee100f5e8dc675a

4
matlab/posterior_sampler_core.m
View File

@ -36,7 +36,7 @@ function myoutput = posterior_sampler_core(myinputs,fblck,nblck,whoiam, ThisMatl
% See the comments in the posterior_sampler.m funtion.
% Copyright © 2006-2017 Dynare Team
% Copyright © 2006-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -86,8 +86,6 @@ oo_ = myinputs.oo_;
if whoiam
% initialize persistent variables in priordens()
priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, bayestopt_.p3,bayestopt_.p4,1);
% initialize persistent variables in prior_draw()
prior_draw(bayestopt_,options_.prior_trunc);
end
MetropolisFolder = CheckPath('metropolis',M_.dname);

13
matlab/posterior_sampler_initialization.m
View File

@ -38,7 +38,7 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npa
% SPECIAL REQUIREMENTS
% None.
% Copyright © 2006-2017 Dynare Team
% Copyright © 2006-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -113,7 +113,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
fprintf(fidlog,[' Number of simulations per block: ' int2str(nruns(1)) '\n']);
fprintf(fidlog,' \n');
if isempty(d)
prior_draw(bayestopt_,options_.prior_trunc);
Prior = dprior(bayestopt_, options_.prior_trunc);
end
if options_.mh_initialize_from_previous_mcmc.status
PreviousFolder0 = options_.mh_initialize_from_previous_mcmc.directory;
@ -132,7 +132,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'];
@ -195,7 +195,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
trial = 1;
while validate == 0 && trial <= 10
if isempty(d)
candidate = prior_draw();
candidate = Prior.draw();
else
if isfield(options_,'mh_init_scale')
if trial==1
@ -228,7 +228,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 +301,6 @@ 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;

17
matlab/prior_posterior_statistics_core.m
View File

@ -30,7 +30,7 @@ function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMa
% SPECIAL REQUIREMENTS.
% None.
% Copyright © 2005-2020 Dynare Team
% Copyright © 2005-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -109,8 +109,8 @@ if ~strcmpi(type,'prior')
logpost=myinputs.logpost;
end
else
x=(prior_draw(bayestopt_,options_.prior_trunc))';
options_=select_qz_criterium_value(options_);
Prior = dprior(bayestopt_,options_.prior_trunc);
options_ = select_qz_criterium_value(options_);
end
if whoiam
Parallel=myinputs.Parallel;
@ -180,9 +180,9 @@ end
if naK
stock_filter_step_ahead =NaN(length(options_.filter_step_ahead),endo_nbr,gend+max(options_.filter_step_ahead),MAX_naK);
end
stock_param = NaN(MAX_nruns,size(x,2));
stock_logpo = NaN(MAX_nruns,1);
stock_ys = NaN(MAX_nruns,endo_nbr);
stock_param = NaN(MAX_nruns, length(bayestopt_.p6));
stock_logpo = NaN(MAX_nruns, 1);
stock_ys = NaN(MAX_nruns, endo_nbr);
if filter_covariance
stock_filter_covariance = zeros(endo_nbr,endo_nbr,gend+1,MAX_filter_covariance);
end
@ -196,7 +196,8 @@ for b=fpar:B
iter=1;
logpo=[];
while (isempty(logpo) || isinf(logpo)) && iter<1000
[deep, logpo] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
deep = Prior.draw();
logpo = evaluate_posterior_kernel(deep, M_, estim_params_, oo_, options_, bayestopt_)
iter=iter+1;
end
if iter==1000
@ -214,7 +215,7 @@ for b=fpar:B
if run_smoother
[dr,info,M_,oo_] =compute_decision_rules(M_,opts_local,oo_);
if ismember(info(1),[3,4])
if ismember(info(1),[3,4])
opts_local.qz_criterium = 1 + (opts_local.qz_criterium-1)*10; %loosen tolerance, useful for big models where BK conditions were tested on restricted state space
[dr,info,M_,oo_] =compute_decision_rules(M_,opts_local,oo_);
end

14
matlab/prior_sampler.m
View File

@ -13,7 +13,7 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_
% SPECIAL REQUIREMENTS
% none
% Copyright © 2009-2018 Dynare Team
% Copyright © 2009-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -31,7 +31,7 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% Initialization.
prior_draw(bayestopt_, options_.prior_trunc);
Prior = dprior(bayestopt_, options_.prior_trunc);
PriorDirectoryName = CheckPath('prior/draws',M_.dname);
work = ~drsave;
iteration = 0;
@ -94,10 +94,10 @@ while iteration < NumberOfSimulations
dyn_waitbar(iteration/NumberOfSimulations,hh,'Please wait. Prior sampler...');
end
loop_indx = loop_indx+1;
params = prior_draw();
M_ = set_all_parameters(params,estim_params_,M_);
[T,R,~,INFO,M_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
dr=oo_.dr;
params = Prior.draw();
M_ = set_all_parameters(params, estim_params_, M_);
[T, R, ~, INFO, M_, oo_] = dynare_resolve(M_, options_, oo_, 'restrict');
dr = oo_.dr;
if ~INFO(1)
INFO=endogenous_prior_restrictions(T,R,M_,options_,oo_);
end
@ -200,4 +200,4 @@ m1 = m0 + (newobs'-m0)/iter;
qq = m1*m1';
s1 = s0 + ( (newobs'*newobs-qq-s0) + (iter-1)*(m0*m0'-qq') )/iter;
mu = m1;
sigma = s1;
sigma = s1;

8
matlab/slice_sampler.m
View File

@ -24,7 +24,7 @@ function [theta, fxsim, neval] = slice_sampler(objective_function,theta,thetapri
% SPECIAL REQUIREMENTS
% none
% Copyright © 2015-2017 Dynare Team
% Copyright © 2015-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -58,6 +58,8 @@ neval = zeros(npar,1);
% % % fname0=fname;
fname = [ int2str(sampler_options.curr_block)];
Prior = dprior(varargin{6},varargin{3}.prior_trunc);
it=0;
while it<npar
it=it+1;
@ -80,8 +82,8 @@ while it<npar
icount=0;
while ~isfinite(fxold) && icount<1000
icount=icount+1;
theta = prior_draw();
if all(theta(:) >= thetaprior(:,1)) && all(theta(:) <= thetaprior(:,2))
theta = Prior.draw();
if all(theta >= thetaprior(:,1)) && all(theta <= thetaprior(:,2))
fxold = -feval(objective_function,theta,varargin{:});
end
end