Commit 0bb00c71 authored by stepan's avatar stepan
Browse files

v4.1::

+ Added new file for prior sampling (prior_sample.m). 
+ Undo previous change on prior_draw.m (transposition of pdraw).
+ Bug correction in prior_draw.m.
+ Added field in options_: the default number of mc simulations in the
prior distribution.


git-svn-id: https://www.dynare.org/svn/dynare/trunk@2431 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 5f02e69a
function get_prior_info(info)
% function dynare_estimation_1(var_list_,dname)
% runs the estimation of the model
%
% INPUTS
% none
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2009 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_
if ~nargin
info = 0;
end
[xparam1,estim_params_,bayestopt_,lb,ub,M_] = set_prior(estim_params_,M_,options_);
plot_priors(bayestopt_,M_,options_);
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+');
% 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_.prior_trunc = prior_trunc_backup ;
for i=1:size(bayestopt_.name,1)
[tmp,TexName] = get_the_name(i,1);
PriorShape = PriorNames{ bayestopt_.pshape(i) };
PriorMean = bayestopt_.pmean(i);
PriorStandardDeviation = bayestopt_.pstdev(i);
switch bayestopt_.pshape(i)
case { 1 , 5 }
LowerBound = bayestopt_.p3(i);
UpperBound = bayestopt_.p4(i);
case { 2 , 4 , 6 }
LowerBound = bayestopt_.p3(i);
UpperBound = '$\infty$';
case 3
if isinf(bayestopt_.p3(i))
LowerBound = '$-\infty$';
else
LowerBound = bayestopt_.p3(i);
end
if isinf(bayestopt_.p4(i))
UpperBound = '$\infty$';
else
UpperBound = bayestopt_.p4(i);
end
otherwise
error('get_prior_info:: Dynare bug!')
end
format_string = build_format_string(bayestopt_,i);
fprintf(fidTeX,format_string, ...
TexName, ...
PriorShape, ...
PriorMean, ...
PriorStandardDeviation, ...
LowerBound, ...
UpperBound, ...
PriorIntervals(i,1), ...
PriorIntervals(i,2) );
end
fclose(fidTeX);
end
M_.dname = M_.fname;
if info% Prior simulations.
results = prior_sampler(1,M_,bayestopt_,options_,oo_);
results.prior.mass
end
function format_string = build_format_string(bayestopt,i)
format_string = ['%s & %s & %6.4f &'];
if isinf(bayestopt.pstdev(i))
format_string = [ format_string , ' %s &'];
else
format_string = [ format_string , ' %6.4f &'];
end
if isinf(bayestopt.p3(i))
format_string = [ format_string , ' %s &'];
else
format_string = [ format_string , ' %6.4f &'];
end
if isinf(bayestopt.p4(i))
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
......@@ -205,3 +205,5 @@ function global_initialization()
options_.homotopy_mode = 0;
options_.homotopy_steps = 1;
% prior analysis
options_.prior_mc = 20000;
\ No newline at end of file
......@@ -45,13 +45,7 @@ if nargin>0 && init
p4 = prior_structure.p4;
number_of_estimated_parameters = length(p1);
a = NaN(number_of_estimated_parameters,1);
b = NaN(number_of_estimated_parameters,1);
if ~isempty(cc)
bounds = cc;
nobounds_flag = 0;
else
nobounds_flag = 1;
end
b = NaN(number_of_estimated_parameters,1);
beta_index = find(prior_shape==1);
gamma_index = find(prior_shape==2);
gaussian_index = find(prior_shape==3);
......@@ -97,6 +91,4 @@ end
if ~isempty(inverse_gamma_2_index)
pdraw(inverse_gamma_2_index) = ...
1./gamrnd(p2(inverse_gamma_2_index)/2,2./p1(inverse_gamma_2_index));
end
pdraw = transpose(pdraw);
\ No newline at end of file
end
\ No newline at end of file
function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_)
% This function builds a (big) prior sample.
%
% INPUTS
% drsave [integer] Scalar. If equal to 1, then dr structure is saved with each prior draw.
% M_ [structure] Model description.
% bayestopt_ [structure] Prior distribution description.
% options_ [structure] Global options of Dynare.
%
% OUTPUTS:
% results [structure] Various statistics.
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2009 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/>.
% Initialization.
prior_draw(1,bayestopt_);
PriorDirectoryName = CheckPath('prior/draws');
work = ~drsave;
iteration = 1;
loop_indx = 0;
file_indx = [];
count_bk_indeterminacy = 0;
count_bk_unstability = 0;
count_bk_singularity = 0;
count_static_var_def = 0;
count_no_steadystate = 0;
count_dll_problem = 0;
count_unknown_problem = 0;
NumberOfSimulations = options_.prior_mc;
NumberOfParameters = length(bayestopt_.pmean);
NumberOfEndogenousVariables = size(M_.endo_names,1);
NumberOfElementsPerFile = ceil(options_.MaxNumberOfBytes/NumberOfParameters/NumberOfEndogenousVariables/8)
if NumberOfSimulations <= NumberOfElementsPerFile
TableOfInformations = [ 1 , NumberOfSimulations , 1] ;
else
NumberOfFiles = fix(NumberOfSimulations/NumberOfElementsPerFile) ;
NumberOfElementsInTheLastFile = NumberOfSimulations - NumberOfElementsPerFile*NumberOfFiles ;
if ~isint(NumberOfSimulations/NumberOfElementsPerFile)
NumberOfFiles = NumberOfFiles + 1 ;
end
TableOfInformations = NaN(NumberOfFiles,3);
TableOfInformations(:,1) = transpose(1:NumberOfFiles) ;
TableOfInformations(1:NumberOfFiles-1,2) = NumberOfElementsPerFile*ones(NumberOfFiles-1,1) ;
TableOfInformations(NumberOfFiles,2) = NumberOfElementsInTheLastFile ;
TableOfInformations(1,3) = 1;
TableOfInformations(2:end,3) = cumsum(TableOfInformations(2:end,2))+1;
end
pdraws = cell(TableOfInformations(1,2),drsave+1) ;
sampled_prior_expectation = zeros(NumberOfParameters,1);
sampled_prior_covariance = zeros(NumberOfParameters,NumberOfParameters);
% Simulations.
while iteration <= NumberOfSimulations
loop_indx = loop_indx+1;
file_indx = find(TableOfInformations(:,3)==iteration);
if ~isempty(file_indx) && file_indx>1
save([ PriorDirectoryName '/prior_draws' int2str(file_indx-1) '.mat' ],'pdraws');
pdraws = cell(TableOfInformations(file_indx,2),drsave+1);
end
params = prior_draw();
set_all_parameters(params);
[dr,INFO] = resol(oo_.steady_state,work);
switch INFO(1)
case 0
pdraws(iteration,1) = {params};
if drsave
pdraws(iteration,2) = {dr};
end
iteration = iteration+1;
[sampled_prior_expectation,sampled_prior_covariance] = ...
recursive_prior_moments(sampled_prior_expectation,sampled_prior_covariance,params,iteration) ;
case 1
count_static_undefined = count_static_undefined + 1;
case 2
count_dll_problem = count_dll_problem + 1;
case 3
count_bk_unstability = count_bk_unstability + 1 ;
case 4
count_bk_indeterminacy = count_bk_indeterminacy + 1 ;
case 5
count_bk_singularity = count_bk_singularity + 1 ;
case 20
count_no_steadystate = count_no_steadystate + 1 ;
otherwise
% To be checked...
count_unknown_problem = count_unknown_problem + 1 ;
continue
end
end
% Save last prior_draw*.mat file
save([ PriorDirectoryName '/prior_draws' int2str(TableOfInformations(end,1)) '.mat' ],'pdraws');
% Get informations about BK conditions and other things...
results.bk.indeterminacy_share = count_bk_indeterminacy/loop_indx;
results.bk.unstability_share = count_bk_unstability/loop_indx;
results.bk.singularity_share = count_bk_singularity/loop_indx;
results.dll.problem_share = count_dll_problem/loop_indx;
results.ss.problem_share = count_no_steadystate/loop_indx;
results.garbage_share = results.bk.indeterminacy_share + ...
results.bk.unstability_share + ...
results.bk.singularity_share + ...
results.dll.problem_share + ...
results.ss.problem_share + ...
(count_unknown_problem/loop_indx) ;
results.prior.mean = sampled_prior_expectation;
results.prior.variance = sampled_prior_covariance;
results.prior.mass = 1-results.garbage_share;
function [mu,sigma] = recursive_prior_moments(m0,s0,newobs,iter)
% Recursive estimation of order one and two moments (expectation and
% covariance matrix). newobs should be a row vector. I do not use the
% function recursive_moments here, because this function is to be used when
% newobs is a 2D array.
m1 = m0 + (newobs'-m0)/iter;
qq = m1*m1';
s1 = s0 + ( (newobs'*newobs-qq-s0) + (iter-1)*(m0*m0'-qq') )/iter;
mu = m1;
sigma = s1;
\ No newline at end of file
function pstat = prior_statistics(info,M_,bayestopt_,options_)
% This function computes various statistics associated with the prior distribution.
%
% INPUTS
% info [integer] 1*p vector.
% M_ [structure] Model description.
% bayestopt_ [structure] Prior distribution description.
% options_ [structure] Global options of Dynare.
%
% OUTPUTS:
% pstat [structure] Various statistics
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2009 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/>.
if ~info(0)
results = prior_sampler(1,M_,bayestopt_,options_);
results.prior.mass
end
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment