+ Various bug fixes related to prior sampling.

+ Removed globals in set_stationary_variables_list.m.


git-svn-id: https://www.dynare.org/svn/dynare/trunk@2771 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
stepan 2009-06-15 14:36:30 +00:00
parent 2d314b32b9
commit f039875f60
13 changed files with 147 additions and 152 deletions

View File

@ -26,9 +26,10 @@ function [info,description] = check_prior_analysis_data(type,M_)
disp('check_prior_analysis_data:: Can''t find any prior draws file!')
return
end
prior_draws_info = dir([ M_.dname '/prior/draws/prior_draws*.mat']);
name_of_the_last_prior_draw_file = prior_draws_info(end).name ;%mhname
date_of_the_last_prior_draw_file = prior_draws_info(end).datenum ;%mhdate
name_of_the_last_prior_draw_file = prior_draws_info(end).name;
date_of_the_last_prior_draw_file = prior_draws_info(end).datenum;
%% Get informations about _posterior_draws files.
if isempty(prior_draws_info)
@ -65,9 +66,9 @@ function [info,description] = check_prior_analysis_data(type,M_)
case 'conditional decomposition'
generic_prior_data_file_name = 'PriorConditionalVarianceDecomposition';
otherwise
disp(['This feature is not yet implemented!')
disp(['This feature is not yet implemented!'])
end
CheckPath('prior/moments')
CheckPath('prior/moments');
pdfinfo = dir([ M_.dname '/prior/' generic_prior_data_file_name '*']);
if isempty(pdfinfo)
info = 4;

View File

@ -4,6 +4,7 @@ function oo_ = compute_moments_varendo(type,options_,M_,oo_,var_list_)
% var_list_. The results are saved in oo_
%
% INPUTS:
% type [string] 'posterior' or 'prior'
% options_ [structure] Dynare structure.
% M_ [structure] Dynare structure (related to model definition).
% oo_ [structure] Dynare structure (results).
@ -34,8 +35,17 @@ function oo_ = compute_moments_varendo(type,options_,M_,oo_,var_list_)
if strcmpi(type,'posterior')
posterior = 1;
if nargin==4
var_list_ = options_.varobs;
end
elseif strcmpi(type,'prior')
posterior = 0;
if nargin==4
var_list_ = options_.prior_analysis_endo_var_list;
if isempty(var_list_)
options_.prior_analysis_var_list = options_.varobs;
end
end
else
disp('compute_moments_varendo:: Unknown type!')
error()

View File

@ -88,11 +88,11 @@ function oo_ = covariance_mc_analysis(NumberOfSimulations,type,dname,fname,varta
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.deciles.' name ' = p_deciles;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.density.' name ' = density;']);
else
eval(['oo_.' NAME 'TheoreticalMoments.dsge.covariance.mean.' name ' = NaN;']);
eval(['oo_.' NAME 'TheoreticalMoments.dsge.covariance.median.' name ' = NaN;']);
eval(['oo_.' NAME 'TheoreticalMoments.dsge.covariance.variance.' name ' = NaN;']);
eval(['oo_.' NAME 'TheoreticalMoments.dsge.covariance.hpdinf.' name ' = NaN;']);
eval(['oo_.' NAME 'TheoreticalMoments.dsge.covariance.hpdsup.' name ' = NaN;']);
eval(['oo_.' NAME 'TheoreticalMoments.dsge.covariance.deciles.' name ' = NaN;']);
eval(['oo_.' NAME 'TheoreticalMoments.dsge.covariance.density.' name ' = NaN;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.mean.' name ' = NaN;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.median.' name ' = NaN;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.variance.' name ' = NaN;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.hpdinf.' name ' = NaN;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.hpdsup.' name ' = NaN;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.deciles.' name ' = NaN;']);
eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.density.' name ' = NaN;']);
end

View File

@ -33,13 +33,6 @@ function [nvar,vartan,NumberOfConditionalDecompFiles] = ...
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Set varlist (vartan)
[ivar,vartan] = set_stationary_variables_list;
nvar = length(ivar);
% Set the size of the auto-correlation function to zero.
nar = options_.ar;
options_.ar = 0;
% Get informations about the _posterior_draws files.
if strcmpi(type,'posterior')
@ -53,6 +46,26 @@ else
error()
end
% Set varlist (vartan)
if ~posterior
if isfield(options_,'varlist')
temp = options_.varlist;
end
options_.varlist = options_.prior_analysis_endo_var_list;
end
[ivar,vartan, options_] = set_stationary_variables_list(options_, M_);
if ~posterior
if exist('temp','var')
options_.varlist = temp;
end
end
nvar = length(ivar);
% Set the size of the auto-correlation function to zero.
nar = options_.ar;
options_.ar = 0;
NumberOfDrawsFiles = length(DrawsFiles);
NumberOfDrawsFiles = rows(DrawsFiles);
@ -65,7 +78,7 @@ if SampleSize<=MaXNumberOfConditionalDecompLines
else
Conditional_decomposition_array = zeros(nvar*(nvar+1)/2,length(Steps),M_.exo_nbr,MaXNumberOfConditionalDecompLines);
NumberOfLinesInTheLastConditionalDecompFile = mod(SampleSize,MaXNumberOfConditionalDecompLines);
NumberOfConditionalDecompFiles = ceil(SampleSize/MaXNumberOfCOnditionalDecompLines);
NumberOfConditionalDecompFiles = ceil(SampleSize/MaXNumberOfConditionalDecompLines);
end
NumberOfConditionalDecompLines = rows(Conditional_decomposition_array);

View File

@ -33,14 +33,6 @@ function [nvar,vartan,CorrFileNumber] = dsge_simulated_theoretical_correlation(S
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
nodecomposition = 1;
% Set varlist (vartan)
[ivar,vartan] = set_stationary_variables_list;
nvar = length(ivar);
% Set the size of the auto-correlation function to nar.
oldnar = options_.ar;
options_.ar = nar;
% Get informations about the _posterior_draws files.
if strcmpi(type,'posterior')
@ -54,6 +46,25 @@ else
error()
end
NumberOfDrawsFiles = length(DrawsFiles);
% Set varlist (vartan)
if ~posterior
if isfield(options_,'varlist')
temp = options_.varlist;
end
options_.varlist = options_.prior_analysis_endo_var_list;
end
[ivar,vartan, options_] = set_stationary_variables_list(options_, M_);
if ~posterior
if exist('temp','var')
options_.varlist = temp;
end
end
nvar = length(ivar);
% Set the size of the auto-correlation function to nar.
oldnar = options_.ar;
options_.ar = nar;
% Number of lines in posterior data files.
MaXNumberOfCorrLines = ceil(options_.MaxNumberOfBytes/(nvar*nvar*nar)/8);

View File

@ -34,14 +34,6 @@ function [nvar,vartan,CovarFileNumber] = dsge_simulated_theoretical_covariance(S
nodecomposition = 1;
% Set varlist (vartan)
[ivar,vartan] = set_stationary_variables_list;
nvar = length(ivar);
% Set the size of the auto-correlation function to zero.
nar = options_.ar;
options_.ar = 0;
% Get informations about the _posterior_draws files.
if strcmpi(type,'posterior')
DrawsFiles = dir([M_.dname '/metropolis/' M_.fname '_' type '_draws*' ]);
@ -55,6 +47,25 @@ else
end
NumberOfDrawsFiles = length(DrawsFiles);
% Set varlist (vartan)
if ~posterior
if isfield(options_,'varlist')
temp = options_.varlist;
end
options_.varlist = options_.prior_analysis_endo_var_list;
end
[ivar,vartan] = set_stationary_variables_list(options_,M_);
if ~posterior
if exist('temp','var')
options_.varlist = temp;
end
end
nvar = length(ivar);
% Set the size of the auto-correlation function to zero.
nar = options_.ar;
options_.ar = 0;
% Number of lines in posterior data files.
MaXNumberOfCovarLines = ceil(options_.MaxNumberOfBytes/(nvar*(nvar+1)/2)/8);

View File

@ -35,14 +35,6 @@ function [nvar,vartan,NumberOfDecompFiles] = ...
nodecomposition = 0;
% Set varlist (vartan)
[ivar,vartan] = set_stationary_variables_list;
nvar = length(ivar);
% Set the size of the auto-correlation function to zero.
nar = options_.ar;
options_.ar = 0;
% Get informations about the _posterior_draws files.
if strcmpi(type,'posterior')
DrawsFiles = dir([M_.dname '/metropolis/' M_.fname '_' type '_draws*' ]);
@ -56,6 +48,27 @@ else
end
NumberOfDrawsFiles = length(DrawsFiles);
% Set varlist (vartan)
if ~posterior
if isfield(options_,'varlist')
temp = options_.varlist;
end
options_.varlist = options_.prior_analysis_endo_var_list;
end
[ivar,vartan,options_] = set_stationary_variables_list(options_,M_);
if ~posterior
if exist('temp','var')
options_.varlist = temp;
end
end
nvar = length(ivar);
% Set the size of the auto-correlation function to zero.
nar = options_.ar;
options_.ar = 0;
nexo = M_.exo_nbr;
NumberOfDrawsFiles = rows(DrawsFiles);

View File

@ -1,9 +1,8 @@
function get_prior_info(info)
% function dynare_estimation_1(var_list_,dname)
% runs the estimation of the model
% Computes various prior statistics.
%
% INPUTS
% none
% info [integer] scalar specifying what has to be done.
%
% OUTPUTS
% none
@ -32,6 +31,8 @@ function get_prior_info(info)
if ~nargin
info = 0;
end
M_.dname = M_.fname;
[xparam1,estim_params_,bayestopt_,lb,ub,M_] = set_prior(estim_params_,M_,options_);
plot_priors(bayestopt_,M_,options_);
@ -94,17 +95,17 @@ function get_prior_info(info)
M_.dname = M_.fname;
if info==1% Prior simulations (BK).
results = prior_sampler(0,M_,bayestopt_,options_,oo_);
% 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)])
results = prior_sampler(0,M_,bayestopt_,options_,oo_);
% 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.
@ -115,7 +116,7 @@ function get_prior_info(info)
look_for_admissible_initial_condition = 1;
scale = 1.0;
iter = 0;
While look_for_admissible_initial_condition
while look_for_admissible_initial_condition
xinit = xparam1+scale*randn(size(xparam1));
if all(xinit>bayestopt_.p3) && all(xinit<bayestopt_.p4)
look_for_admissible_initial_condition = 0;
@ -130,11 +131,11 @@ function get_prior_info(info)
end
% Maximization
[xparams,lpd,hessian] = ...
maximize_prior_density(xinit, bayestopt_.pshape, ...
bayestopt_.p6, ...
bayestopt_.p7, ...
bayestopt_.p3, ...
bayestopt_.p4);
maximize_prior_density(xinit, bayestopt_.pshape, ...
bayestopt_.p6, ...
bayestopt_.p7, ...
bayestopt_.p3, ...
bayestopt_.p4);
% Display the results.
disp(' ')
disp(' ')
@ -142,85 +143,19 @@ function get_prior_info(info)
disp('PRIOR OPTIMIZATION')
disp('------------------')
disp(' ')
for i = 1:length(xparams)
disp(['deep parameter ' int2str(i) ': ' get_the_name(i,0) '.'])
disp([' Initial condition ....... ' num2str(xinit(i)) '.'])
disp([' Prior mode .............. ' num2str(bayestopt_.p5(i)) '.'])
disp([' Optimized prior mode .... ' num2str(xparams(i)) '.'])
disp(' ')
end
for i = 1:length(xparams)
disp(['deep parameter ' int2str(i) ': ' get_the_name(i,0) '.'])
disp([' Initial condition ....... ' num2str(xinit(i)) '.'])
disp([' Prior mode .............. ' num2str(bayestopt_.p5(i)) '.'])
disp([' Optimized prior mode .... ' num2str(xparams(i)) '.'])
disp(' ')
end
end
if info==3% Prior simulations (BK + 2nd order moments).
results = prior_sampler(1,M_,bayestopt_,options_,oo_);
% 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)])
assignin('base','prior_mass_analysis',results);
% Compute prior moments.
PriorMomentsDirectoryName = CheckPath('prior/moments');
prior_draws_info = dir([ M_.dname '/prior/draws/prior_draws*.mat']);
number_of_prior_draws_files = length(prior_draws_info);
total_number_of_simulations = 0;
nvar = rows(options_.prior_analysis_endo_var_list);
if nvar == 0
nvar = M_.endo_nbr;
ivar = [1:nvar]';
else
ivar=zeros(nvar,1);
for i=1:nvar
i_tmp = strmatch(options_.prior_analysis_endo_var_list(i,:),M_.endo_names,'exact');
if isempty(i_tmp)
error (['One of the variable specified does not exist']) ;
else
ivar(i) = i_tmp;
end
end
end
for f=1:number_of_prior_draws_files
load([ M_.dname '/prior/draws/prior_draws' int2str(f) '.mat']);
number_of_simulations = length(pdraws);
total_number_of_simulations = total_number_of_simulations + number_of_simulations;
covariance_cell = cell(number_of_simulations,1);
correlation_cell = cell(number_of_simulations,1);
decomposition_cell = cell(number_of_simulations,1);
for s=1:number_of_simulations
[gamma_y,ivar] = th_autocovariances(pdraws{s,2},ivar,M_,options_);
covariance_cell(s) = {vech(gamma_y{1})};
tmp = zeros(length(ivar),options_.ar);
for i=1:length(ivar)
for lag=1:options_.ar
tmp(i,lag) = gamma_y{i,lag+1};
end
end
correlation_cell(s) = {tmp};
decomposition_cell(s) = {gamma_y{options_.ar+2}};
end
save([ PriorMomentsDirectoryName '/prior_moments_draws' int2str(f) '.mat' ],'covariance_cell','correlation_cell','decomposition_cell');
end
clear('covariance_cell','correlation_cell','decomposition_cell')
prior_moments_info = dir([ M_.dname '/prior/moments/prior_moments*.mat']);
number_of_prior_moments_files = length(prior_moments_info);
% Covariance analysis
disp(' ')
disp('-------------------------')
disp('Prior variance analysis')
disp('-------------------------')
disp(' ')
for i=1:length(ivar)
for file = 1:number_of_prior_moments_file
load()
end
end
if info==3% Prior simulations (2nd order moments).
oo_ = compute_moments_varendo('prior',options_,M_,oo_);
end
function format_string = build_format_string(bayestopt,i)
format_string = ['%s & %s & %6.4f &'];

View File

@ -217,4 +217,5 @@ function global_initialization()
options_.homotopy_steps = 1;
% prior analysis
options_.prior_mc = 20000;
options_.prior_mc = 20000;
options_.prior_analysis_endo_var_list = [];

View File

@ -33,7 +33,7 @@ function oo_ = posterior_analysis(type,arg1,arg2,arg3,options_,M_,oo_)
case {4,5}
oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
case 6
[ivar,vartan] = set_stationary_variables_list;
[ivar,vartan] = set_stationary_variables_list(options_,M_);
nvar = length(ivar);
oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_,nvar,vartan);
otherwise

View File

@ -34,7 +34,7 @@ function oo_ = prior_analysis(type,arg1,arg2,arg3,options_,M_,oo_)
case {4,5}
oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_);
case 6
[ivar,vartan] = set_stationary_variables_list;
[ivar,vartan] = set_stationary_variables_list(options_,M_);
nvar = length(ivar);
oo_ = job(type,SampleSize,arg1,arg2,arg3,options_,M_,oo_,nvar,vartan);
otherwise

View File

@ -1,22 +1,23 @@
function [ivar,vartan] = set_stationary_variables_list()
function [ivar,vartan,options_] = set_stationary_variables_list(options_,M_)
% This function builds a vector of indices targeting to the stationary
% variables in varlist.
%
% INPUTS
% None.
%
% o options_ [structure] Describes global options.
% o M_ [structure] Describes the model.
% OUTPUTS
% o ivar [integer] nvar*1 vector of indices (nvar is the number
% of stationary variables).
% o vartan [char] array of characters (with nvar rows).
%
% o ivar [integer] nvar*1 vector of indices (nvar is the number
% of stationary variables).
% o vartan [char] array of characters (with nvar rows).
% o options_ [structure] Describes global options.
%
% ALGORITHM
% None.
%
% SPECIAL REQUIREMENTS
% None.
% Copyright (C) 2007 Dynare Team
% Copyright (C) 2007-2009 Dynare Team
%
% This file is part of Dynare.
%
@ -33,7 +34,6 @@ function [ivar,vartan] = set_stationary_variables_list()
% 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_
varlist = options_.varlist;
if isempty(varlist)
varlist = options_.varobs;

View File

@ -76,7 +76,7 @@ function oo_ = variance_decomposition_mc_analysis(NumberOfSimulations,type,dname
p_median = t1;
p_var = 0;
hpd_interval = NaN(2,1);
post_deciles = NaN(9,1);
p_deciles = NaN(9,1);
density = NaN;
else
[p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...