Skip to content
Snippets Groups Projects
Commit 73cdaa30 authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files

Replaced get_prior_info routine by prior routine.

This new routine gives the possibility to obtain a description of the
prior from the command line.

>> prior table

will display a table with the prior mean, prior mode, prior std, prior
bounds and prior HPD interval.

>> prior plot

will plot the prior densities.

>> prior optimize

will trigger the maximization of the prior density and display the
results (often, because of BK conditions or other issues, Dynare is even
not able to get the prior mode).

>> prior moments

will display the moments of the endogenous variables at the prior mode.

>> prior simulate

will run a Monte-Carlo, by sampling from the prior and discarding vector
of parameters such that the steady state does not exist or the BK
conditions are not met, and return an estimate of the effective prior
mass.
parent daaed377
Branches
Tags
No related merge requests found
function varargout = prior(varargin)
% Computes various prior statistics and display them in the command window.
%
% INPUTS
% 'table', 'moments', 'optimize', 'simulate', 'plot'
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2015 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 <http://www.gnu.org/licenses/>.
global options_ M_ estim_params_ bayestopt_ oo_ objective_function_penalty_base
donesomething = false;
% Temporarly change qz_criterium option value
changed_qz_criterium_flag = 0;
if isempty(options_.qz_criterium)
options_.qz_criterium = 1+1e-9;
changed_qz_criterium_flag = 1;
end
M_.dname = M_.fname;
% Temporarly set options_.order equal to one
order = options_.order;
options_.order = 1;
[xparam1,estim_params_,bayestopt_,lb,ub,M_] = set_prior(estim_params_,M_,options_);
if ismember('plot', varargin)
plot_priors(bayestopt_,M_,estim_params_,options_)
donesomething = true;
end
if ismember('table', varargin)
print_table_prior(lb, ub, options_, M_, bayestopt_, estim_params_);
donesomething = true;
end
if ismember('simulate', varargin) % Prior simulations (BK).
results = prior_sampler(0,M_,bayestopt_,options_,oo_,estim_params_);
% Display prior mass info
skipline(2)
disp(['Prior mass = ' num2str(results.prior.mass)])
disp(['BK indeterminacy share = ' num2str(results.bk.indeterminacy_share)])
disp(['BK unstability share = ' num2str(results.bk.unstability_share)])
disp(['BK singularity share = ' num2str(results.bk.singularity_share)])
disp(['Complex jacobian share = ' num2str(results.jacobian.problem_share)])
disp(['mjdgges crash share = ' num2str(results.dll.problem_share)])
disp(['Steady state problem share = ' num2str(results.ss.problem_share)])
disp(['Complex steady state share = ' num2str(results.ss.complex_share)])
disp(['Analytical steady state problem share = ' num2str(results.ass.problem_share)])
skipline(2)
donesomething = true;
end
if ismember('optimize', varargin) % Prior optimization.
optimize_prior(options_, M_, oo_, bayestopt_, estim_params_);
donesomething = true;
end
if ismember('moments', varargin) % Prior simulations (2nd order moments).
% Set estimated parameters to the prior mode...
xparam1 = bayestopt_.p5;
% ... Except for uniform priors!
k = find(isnan(xparam1));
xparam1(k) = bayestopt_.p5(k);
% Update vector of parameters and covariance matrices
M_ = set_all_parameters(xparam1,estim_params_,M_);
% Check model.
check_model(M_);
% Compute state space representation of the model.
oo_.dr=set_state_space(oo_.dr, M_, options_);
% Solve model
[dr,info, M_ ,options_ , oo_] = resol(0, M_ , options_ ,oo_);
% Compute and display second order moments
disp_th_moments(oo_.dr,[]);
skipline(2)
donesomething = true;
end
if changed_qz_criterium_flag
options_.qz_criterium = [];
end
options_.order = order;
if ~donesomething
error('prior: Unexpected arguments!')
end
function format_string = build_format_string(PriorStandardDeviation,LowerBound,UpperBound)
format_string = ['%s & %s & %6.4f &'];
if ~isnumeric(PriorStandardDeviation)
format_string = [ format_string , ' %s &'];
else
format_string = [ format_string , ' %6.4f &'];
end
if ~isnumeric(LowerBound)
format_string = [ format_string , ' %s &'];
else
format_string = [ format_string , ' %6.4f &'];
end
if ~isnumeric(UpperBound)
format_string = [ format_string , ' %s &'];
else
format_string = [ format_string , ' %6.4f &'];
end
format_string = [ format_string , ' %6.4f & %6.4f \\\\ \n'];
\ No newline at end of file
...@@ -58,6 +58,7 @@ addpath([dynareroot '/parallel/']) ...@@ -58,6 +58,7 @@ addpath([dynareroot '/parallel/'])
addpath([dynareroot '/particle/']) addpath([dynareroot '/particle/'])
addpath([dynareroot '/gsa/']) addpath([dynareroot '/gsa/'])
addpath([dynareroot '/ep/']) addpath([dynareroot '/ep/'])
addpath([dynareroot '/cli/'])
addpath([dynareroot '/lmmcp/']) addpath([dynareroot '/lmmcp/'])
addpath([dynareroot '/optimization/']) addpath([dynareroot '/optimization/'])
addpath([dynareroot '/modules/dates/src/']) addpath([dynareroot '/modules/dates/src/'])
......
function results = get_prior_info(info,plt_flag)
% Computes various prior statistics.
%
% INPUTS
% info [integer] scalar specifying what has to be done.
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2009-2012 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 <http://www.gnu.org/licenses/>.
global options_ M_ estim_params_ oo_ objective_function_penalty_base
if ~nargin
info = 0;
plt_flag = 0;
end
if nargin==1
plt_flag = 1;
end
% Initialize returned variable.
results = [];
changed_qz_criterium_flag = 0;
if isempty(options_.qz_criterium)
options_.qz_criterium = 1+1e-9;
changed_qz_criterium_flag = 1;
end
M_.dname = M_.fname;
% Temporarly set options_.order equal to one
order = options_.order;
options_.order = 1;
[xparam1,estim_params_,bayestopt_,lb,ub,M_] = set_prior(estim_params_,M_,options_);
if plt_flag
plot_priors(bayestopt_,M_,estim_params_,options_)
end
PriorNames = { 'Beta' , 'Gamma' , 'Gaussian' , 'Inverted Gamma' , 'Uniform' , 'Inverted Gamma -- 2' };
if size(M_.param_names,1)==size(M_.param_names_tex,1)% All the parameters have a TeX name.
fidTeX = fopen('priors_data.tex','w+');
fprintf(fidTeX,'%% TeX-table generated by get_prior_info (Dynare).\n');
fprintf(fidTeX,'%% Prior Information\n');
fprintf(fidTeX,['%% ' datestr(now,0)]);
fprintf(fidTeX,' \n');
fprintf(fidTeX,' \n');
fprintf(fidTeX,'\\begin{center}\n');
fprintf(fidTeX,'\\begin{longtable}{l|ccccccc} \n');
fprintf(fidTeX,'\\caption{Prior information (parameters)}\\\\\n ');
fprintf(fidTeX,'\\label{Table:Prior}\\\\\n');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Lower Bound & Upper Bound & LB Untrunc. 80\\%% HPDI & UB Untrunc. 80\\%% HPDI \\\\ \n');
fprintf(fidTeX,'\\hline \\endfirsthead \n');
fprintf(fidTeX,'\\caption{(continued)}\\\\\n ');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Lower Bound & Upper Bound & LB Untrunc. 80\\%% HPDI & UB Untrunc. 80\\%% HPDI \\\\ \n');
fprintf(fidTeX,'\\hline \\endhead \n');
fprintf(fidTeX,'\\hline \\multicolumn{8}{r}{(Continued on next page)} \\\\ \\hline \\endfoot \n');
fprintf(fidTeX,'\\hline \\hline \\endlastfoot \n');
% Column 1: a string for the name of the prior distribution.
% Column 2: the prior mean.
% Column 3: the prior standard deviation.
% Column 4: the lower bound of the prior density support.
% Column 5: the upper bound of the prior density support.
% Column 6: the lower bound of the interval containing 80% of the prior mass.
% Column 7: the upper bound of the interval containing 80% of the prior mass.
prior_trunc_backup = options_.prior_trunc ;
options_.prior_trunc = (1-options_.prior_interval)/2 ;
PriorIntervals = prior_bounds(bayestopt_,options_) ;
options_.prior_trunc = prior_trunc_backup ;
for i=1:size(bayestopt_.name,1)
[tmp,TexName] = get_the_name(i,1,M_,estim_params_,options_);
PriorShape = PriorNames{ bayestopt_.pshape(i) };
PriorMean = bayestopt_.p1(i);
PriorStandardDeviation = bayestopt_.p2(i);
switch bayestopt_.pshape(i)
case { 1 , 5 }
LowerBound = bayestopt_.p3(i);
UpperBound = bayestopt_.p4(i);
if ~isinf(lb(i))
LowerBound=max(LowerBound,lb(i));
end
if ~isinf(ub(i))
UpperBound=min(UpperBound,ub(i));
end
case { 2 , 4 , 6 }
LowerBound = bayestopt_.p3(i);
if ~isinf(lb(i))
LowerBound=max(LowerBound,lb(i));
end
if ~isinf(ub(i))
UpperBound=ub(i);
else
UpperBound = '$\infty$';
end
case 3
if isinf(bayestopt_.p3(i)) && isinf(lb(i))
LowerBound = '$-\infty$';
else
LowerBound = bayestopt_.p3(i);
if ~isinf(lb(i))
LowerBound=max(LowerBound,lb(i));
end
end
if isinf(bayestopt_.p4(i)) && isinf(ub(i))
UpperBound = '$\infty$';
else
UpperBound = bayestopt_.p4(i);
if ~isinf(ub(i))
UpperBound=min(UpperBound,ub(i));
end
end
otherwise
error('get_prior_info:: Dynare bug!')
end
format_string = build_format_string(PriorStandardDeviation,LowerBound,UpperBound);
fprintf(fidTeX,format_string, ...
TexName, ...
PriorShape, ...
PriorMean, ...
PriorStandardDeviation, ...
LowerBound, ...
UpperBound, ...
PriorIntervals.lb(i), ...
PriorIntervals.ub(i) );
end
fprintf(fidTeX,'\\end{longtable}\n ');
fprintf(fidTeX,'\\end{center}\n');
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
M_.dname = M_.fname;
if info==1% Prior simulations (BK).
results = prior_sampler(0,M_,bayestopt_,options_,oo_,estim_params_);
% Display prior mass info
disp(['Prior mass = ' num2str(results.prior.mass)])
disp(['BK indeterminacy share = ' num2str(results.bk.indeterminacy_share)])
disp(['BK unstability share = ' num2str(results.bk.unstability_share)])
disp(['BK singularity share = ' num2str(results.bk.singularity_share)])
disp(['Complex jacobian share = ' num2str(results.jacobian.problem_share)])
disp(['mjdgges crash share = ' num2str(results.dll.problem_share)])
disp(['Steady state problem share = ' num2str(results.ss.problem_share)])
disp(['Complex steady state share = ' num2str(results.ss.complex_share)])
disp(['Analytical steady state problem share = ' num2str(results.ass.problem_share)])
end
if info==2% Prior optimization.
% Initialize to the prior mode if possible
oo_.dr=set_state_space(oo_.dr,M_,options_);
k = find(~isnan(bayestopt_.p5));
xparam1(k) = bayestopt_.p5(k);
% Pertubation of the initial condition.
look_for_admissible_initial_condition = 1;
scale = 1.0;
iter = 0;
while look_for_admissible_initial_condition
xinit = xparam1+scale*randn(size(xparam1));
if all(xinit(:)>bayestopt_.p3) && all(xinit(:)<bayestopt_.p4)
M_ = set_all_parameters(xinit,estim_params_,M_);
[dr,INFO,M_,options_,oo_] = resol(0,M_,options_,oo_);
if ~INFO(1)
look_for_admissible_initial_condition = 0;
end
else
if iter == 2000;
scale = scale/1.1;
iter = 0;
else
iter = iter+1;
end
end
end
objective_function_penalty_base = minus_logged_prior_density(xinit, bayestopt_.pshape, ...
bayestopt_.p6, ...
bayestopt_.p7, ...
bayestopt_.p3, ...
bayestopt_.p4,options_,M_,estim_params_,oo_);
% Maximization
[xparams,lpd,hessian] = ...
maximize_prior_density(xinit, bayestopt_.pshape, ...
bayestopt_.p6, ...
bayestopt_.p7, ...
bayestopt_.p3, ...
bayestopt_.p4,options_,M_,bayestopt_,estim_params_,oo_);
% Display the results.
skipline(2)
disp('------------------')
disp('PRIOR OPTIMIZATION')
disp('------------------')
skipline()
for i = 1:length(xparams)
disp(['deep parameter ' int2str(i) ': ' get_the_name(i,0,M_,estim_params_,options_) '.'])
disp([' Initial condition ....... ' num2str(xinit(i)) '.'])
disp([' Prior mode .............. ' num2str(bayestopt_.p5(i)) '.'])
disp([' Optimized prior mode .... ' num2str(xparams(i)) '.'])
skipline()
end
end
if info==3% Prior simulations (2nd order moments).
oo_ = compute_moments_varendo('prior',options_,M_,oo_);
end
if changed_qz_criterium_flag
options_.qz_criterium = [];
end
options_.order = order;
function format_string = build_format_string(PriorStandardDeviation,LowerBound,UpperBound)
format_string = ['%s & %s & %6.4f &'];
if ~isnumeric(PriorStandardDeviation)
format_string = [ format_string , ' %s &'];
else
format_string = [ format_string , ' %6.4f &'];
end
if ~isnumeric(LowerBound)
format_string = [ format_string , ' %s &'];
else
format_string = [ format_string , ' %6.4f &'];
end
if ~isnumeric(UpperBound)
format_string = [ format_string , ' %s &'];
else
format_string = [ format_string , ' %6.4f &'];
end
format_string = [ format_string , ' %6.4f & %6.4f \\\\ \n'];
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment