+ Bug fixes related to the computation of prior (and posterior) moments.

+ Cosmetic changes.


git-svn-id: https://www.dynare.org/svn/dynare/trunk@2791 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
stepan 2009-06-26 10:21:30 +00:00
parent 53b0cf6306
commit d6fe54dc18
7 changed files with 53 additions and 47 deletions

View File

@ -35,6 +35,7 @@ function PackedConditionalVarianceDecomposition = conditional_variance_decomposi
ConditionalVariance = zeros(StateSpaceModel.number_of_state_equations,StateSpaceModel.number_of_state_equations); ConditionalVariance = zeros(StateSpaceModel.number_of_state_equations,StateSpaceModel.number_of_state_equations);
ConditionalVariance = repmat(ConditionalVariance,[1 1 length(Steps) StateSpaceModel.number_of_state_innovations]); ConditionalVariance = repmat(ConditionalVariance,[1 1 length(Steps) StateSpaceModel.number_of_state_innovations]);
BB = StateSpaceModel.impulse_matrix*transpose(StateSpaceModel.impulse_matrix); BB = StateSpaceModel.impulse_matrix*transpose(StateSpaceModel.impulse_matrix);
for h = 1:length(Steps) for h = 1:length(Steps)
for t = 0:Steps(h) for t = 0:Steps(h)
for i=1:StateSpaceModel.number_of_state_innovations for i=1:StateSpaceModel.number_of_state_innovations
@ -44,6 +45,7 @@ function PackedConditionalVarianceDecomposition = conditional_variance_decomposi
end end
end end
end end
ConditionalVariance = ConditionalVariance(SubsetOfVariables,SubsetOfVariables,:,:); ConditionalVariance = ConditionalVariance(SubsetOfVariables,SubsetOfVariables,:,:);
NumberOfVariables = length(SubsetOfVariables); NumberOfVariables = length(SubsetOfVariables);
PackedConditionalVarianceDecomposition = zeros(NumberOfVariables*(NumberOfVariables+1)/2,length(Steps),StateSpaceModel.number_of_state_innovations); PackedConditionalVarianceDecomposition = zeros(NumberOfVariables*(NumberOfVariables+1)/2,length(Steps),StateSpaceModel.number_of_state_innovations);

View File

@ -45,7 +45,7 @@ function oo_ = conditional_variance_decomposition_mc_analysis(NumberOfSimulation
if isfield(temporary_structure,'dsge') if isfield(temporary_structure,'dsge')
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge;']) eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge;'])
if isfield(temporary_structure,'ConditionalVarianceDecomposition') if isfield(temporary_structure,'ConditionalVarianceDecomposition')
eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge.VarianceDecomposition.mean;']) eval(['temporary_structure = oo_.' TYPE 'TheoreticalMoments.dsge.ConditionalVarianceDecomposition.mean;'])
if isfield(temporary_structure,name) if isfield(temporary_structure,name)
if sum(Steps-temporary_structure.(name)(1,:)) == 0 if sum(Steps-temporary_structure.(name)(1,:)) == 0
% Nothing (new) to do here... % Nothing (new) to do here...

View File

@ -40,6 +40,7 @@ if strcmpi(type,'posterior')
posterior = 1; posterior = 1;
elseif strcmpi(type,'prior') elseif strcmpi(type,'prior')
DrawsFiles = dir([M_.dname '/prior/draws/' type '_draws*' ]); DrawsFiles = dir([M_.dname '/prior/draws/' type '_draws*' ]);
CheckPath('prior/moments');
posterior = 0; posterior = 0;
else else
disp('dsge_simulated_theoretical_conditional_variance_decomposition:: Unknown type!') disp('dsge_simulated_theoretical_conditional_variance_decomposition:: Unknown type!')
@ -53,7 +54,7 @@ if ~posterior
end end
options_.varlist = options_.prior_analysis_endo_var_list; options_.varlist = options_.prior_analysis_endo_var_list;
end end
[ivar,vartan, options_] = set_stationary_variables_list(options_, M_); [ivar,vartan ] = set_stationary_variables_list(options_, M_);
if ~posterior if ~posterior
if exist('temp','var') if exist('temp','var')
options_.varlist = temp; options_.varlist = temp;
@ -65,9 +66,6 @@ nvar = length(ivar);
nar = options_.ar; nar = options_.ar;
options_.ar = 0; options_.ar = 0;
NumberOfDrawsFiles = length(DrawsFiles);
NumberOfDrawsFiles = rows(DrawsFiles); NumberOfDrawsFiles = rows(DrawsFiles);
NumberOfSavedElementsPerSimulation = nvar*(nvar+1)/2*M_.exo_nbr*length(Steps); NumberOfSavedElementsPerSimulation = nvar*(nvar+1)/2*M_.exo_nbr*length(Steps);
MaXNumberOfConditionalDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation/8); MaXNumberOfConditionalDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation/8);
@ -81,21 +79,13 @@ else
NumberOfConditionalDecompFiles = ceil(SampleSize/MaXNumberOfConditionalDecompLines); NumberOfConditionalDecompFiles = ceil(SampleSize/MaXNumberOfConditionalDecompLines);
end end
NumberOfConditionalDecompLines = rows(Conditional_decomposition_array); NumberOfConditionalDecompLines = size(Conditional_decomposition_array,4);
ConditionalDecompFileNumber = 1; ConditionalDecompFileNumber = 0;
StateSpaceModel.number_of_state_equations = M_.endo_nbr; StateSpaceModel.number_of_state_equations = M_.endo_nbr;
StateSpaceModel.number_of_state_innovations = M_.exo_nbr; StateSpaceModel.number_of_state_innovations = M_.exo_nbr;
endo_nbr = M_.endo_nbr; first_call = 1;
nstatic = oo_.dr.nstatic;
npred = oo_.dr.npred;
iv = (1:endo_nbr)';
ic = [ nstatic+(1:npred) endo_nbr+(1:size(oo_.dr.ghx,2)-npred) ]';
aux = oo_.dr.transition_auxiliary_variables;
k = find(aux(:,2) > npred);
aux(:,2) = aux(:,2) + nstatic;
aux(k,2) = aux(k,2) + oo_.dr.nfwrd;
linea = 0; linea = 0;
for file = 1:NumberOfDrawsFiles for file = 1:NumberOfDrawsFiles
@ -115,11 +105,28 @@ for file = 1:NumberOfDrawsFiles
set_parameters(pdraws{linee,1}); set_parameters(pdraws{linee,1});
[dr,info] = resol(oo_.steady_state,0); [dr,info] = resol(oo_.steady_state,0);
end end
if first_call
endo_nbr = M_.endo_nbr;
nstatic = dr.nstatic;
npred = dr.npred;
iv = (1:endo_nbr)';
ic = [ nstatic+(1:npred) endo_nbr+(1:size(dr.ghx,2)-npred) ]';
aux = dr.transition_auxiliary_variables;
k = find(aux(:,2) > npred);
aux(:,2) = aux(:,2) + nstatic;
aux(k,2) = aux(k,2) + dr.nfwrd;
StateSpaceModel.number_of_state_equations = M_.endo_nbr+rows(aux);
StateSpaceModel.number_of_state_innovations = M_.exo_nbr;
first_call = 0;
clear('endo_nbr','nstatic','npred','k');
end
[StateSpaceModel.transition_matrix,StateSpaceModel.impulse_matrix] = kalman_transition_matrix(dr,iv,ic,aux,M_.exo_nbr); [StateSpaceModel.transition_matrix,StateSpaceModel.impulse_matrix] = kalman_transition_matrix(dr,iv,ic,aux,M_.exo_nbr);
StateSpaceModel.state_innovations_covariance_matrix = M_.Sigma_e; StateSpaceModel.state_innovations_covariance_matrix = M_.Sigma_e;
clear('dr'); clear('dr');
Conditional_decomposition_array(:,:,:,linea) = conditional_variance_decomposition(StateSpaceModel, Steps, ivar); Conditional_decomposition_array(:,:,:,linea) = conditional_variance_decomposition(StateSpaceModel, Steps, ivar);
if linea == NumberOfConditionalDecompLines if linea == NumberOfConditionalDecompLines
ConditionalDecompFileNumber = ConditionalDecompFileNumber + 1;
linea = 0;
if posterior if posterior
save([M_.dname '/metropolis/' M_.fname '_PosteriorConditionalVarianceDecomposition' int2str(ConditionalDecompFileNumber) '.mat' ], ... save([M_.dname '/metropolis/' M_.fname '_PosteriorConditionalVarianceDecomposition' int2str(ConditionalDecompFileNumber) '.mat' ], ...
'Conditional_decomposition_array'); 'Conditional_decomposition_array');
@ -127,14 +134,10 @@ for file = 1:NumberOfDrawsFiles
save([M_.dname '/prior/moments/' M_.fname '_PriorConditionalVarianceDecomposition' int2str(ConditionalDecompFileNumber) '.mat' ], ... save([M_.dname '/prior/moments/' M_.fname '_PriorConditionalVarianceDecomposition' int2str(ConditionalDecompFileNumber) '.mat' ], ...
'Conditional_decomposition_array'); 'Conditional_decomposition_array');
end end
ConditionalDecompFileNumber = ConditionalDecompFileNumber + 1; if (ConditionalDecompFileNumber==NumberOfConditionalDecompFiles-1)% Prepare last round.
linea = 0; Conditional_decomposition_array = zeros(nvar*(nvar+1)/2, length(Steps),M_.exo_nbr,NumberOfLinesInTheLastConditionalDecompFile) ;
test = ConditionalDecompFileNumber-NumberOfConditionalDecompFiles;
if ~test% Prepare the last round...
Conditional_decomposition_array = zeros(nvar*(nvar+1)/2,length(Steps),M_.exo_nbr,NumberOfLinesInTheLastConditionalDecompFile);
NumberOfConditionalDecompLines = NumberOfLinesInTheLastConditionalDecompFile; NumberOfConditionalDecompLines = NumberOfLinesInTheLastConditionalDecompFile;
ConditionalDecompFileNumber = ConditionalDecompFileNumber - 1; elseif ConditionalDecompFileNumber<NumberOfConditionalDecompFiles-1
elseif test<0;
Conditional_decomposition_array = zeros(nvar*(nvar+1)/2,length(Steps),M_.exo_nbr,MaXNumberOfConditionalDecompLines); Conditional_decomposition_array = zeros(nvar*(nvar+1)/2,length(Steps),M_.exo_nbr,MaXNumberOfConditionalDecompLines);
else else
clear('Conditional_decomposition_array'); clear('Conditional_decomposition_array');

View File

@ -40,6 +40,7 @@ if strcmpi(type,'posterior')
posterior = 1; posterior = 1;
elseif strcmpi(type,'prior') elseif strcmpi(type,'prior')
DrawsFiles = dir([M_.dname '/prior/draws/' type '_draws*' ]); DrawsFiles = dir([M_.dname '/prior/draws/' type '_draws*' ]);
CheckPath('prior/moments');
posterior = 0; posterior = 0;
else else
disp('dsge_simulated_theoretical_correlation:: Unknown type!'); disp('dsge_simulated_theoretical_correlation:: Unknown type!');

View File

@ -40,6 +40,7 @@ if strcmpi(type,'posterior')
posterior = 1; posterior = 1;
elseif strcmpi(type,'prior') elseif strcmpi(type,'prior')
DrawsFiles = dir([M_.dname '/prior/draws/' type '_draws*' ]); DrawsFiles = dir([M_.dname '/prior/draws/' type '_draws*' ]);
CheckPath('prior/moments');
posterior = 0; posterior = 0;
else else
disp('dsge_simulated_theoretical_covariance:: Unknown type!') disp('dsge_simulated_theoretical_covariance:: Unknown type!')

View File

@ -41,6 +41,7 @@ if strcmpi(type,'posterior')
posterior = 1; posterior = 1;
elseif strcmpi(type,'prior') elseif strcmpi(type,'prior')
DrawsFiles = dir([M_.dname '/prior/draws/' type '_draws*' ]); DrawsFiles = dir([M_.dname '/prior/draws/' type '_draws*' ]);
CheckPath('prior/moments');
posterior = 0; posterior = 0;
else else
disp('dsge_simulated_theoretical_variance_decomposition:: Unknown type!') disp('dsge_simulated_theoretical_variance_decomposition:: Unknown type!')

View File

@ -34,7 +34,7 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_)
prior_draw(1,bayestopt_); prior_draw(1,bayestopt_);
PriorDirectoryName = CheckPath('prior/draws'); PriorDirectoryName = CheckPath('prior/draws');
work = ~drsave; work = ~drsave;
iteration = 1; iteration = 0;
loop_indx = 0; loop_indx = 0;
file_indx = []; file_indx = [];
count_bk_indeterminacy = 0; count_bk_indeterminacy = 0;
@ -55,48 +55,39 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_)
if NumberOfSimulations <= NumberOfElementsPerFile if NumberOfSimulations <= NumberOfElementsPerFile
TableOfInformations = [ 1 , NumberOfSimulations , 1] ; TableOfInformations = [ 1 , NumberOfSimulations , 1] ;
else else
NumberOfFiles = fix(NumberOfSimulations/NumberOfElementsPerFile) ; NumberOfFiles = ceil(NumberOfSimulations/NumberOfElementsPerFile) ;
NumberOfElementsInTheLastFile = NumberOfSimulations - NumberOfElementsPerFile*(NumberOfFiles-1) ; NumberOfElementsInTheLastFile = NumberOfSimulations - NumberOfElementsPerFile*(NumberOfFiles-1) ;
if ~isint(NumberOfSimulations/NumberOfElementsPerFile) TableOfInformations = NaN(NumberOfFiles,3) ;
NumberOfFiles = NumberOfFiles + 1 ;
end
TableOfInformations = NaN(NumberOfFiles,3);
TableOfInformations(:,1) = transpose(1:NumberOfFiles) ; TableOfInformations(:,1) = transpose(1:NumberOfFiles) ;
TableOfInformations(1:NumberOfFiles-1,2) = NumberOfElementsPerFile*ones(NumberOfFiles-1,1) ; TableOfInformations(1:NumberOfFiles-1,2) = NumberOfElementsPerFile*ones(NumberOfFiles-1,1) ;
TableOfInformations(NumberOfFiles,2) = NumberOfElementsInTheLastFile ; TableOfInformations(NumberOfFiles,2) = NumberOfElementsInTheLastFile ;
TableOfInformations(1,3) = 1; TableOfInformations(1,3) = 1;
TableOfInformations(2:end,3) = cumsum(TableOfInformations(2:end,2))+1; TableOfInformations(2:end,3) = cumsum(TableOfInformations(1:end-1,2))+1;
end end
pdraws = cell(TableOfInformations(1,2),drsave+1) ; pdraws = cell(TableOfInformations(1,2),drsave+1) ;
sampled_prior_expectation = zeros(NumberOfParameters,1); sampled_prior_expectation = zeros(NumberOfParameters,1);
sampled_prior_covariance = zeros(NumberOfParameters,NumberOfParameters); sampled_prior_covariance = zeros(NumberOfParameters,NumberOfParameters);
file_line_number = 1; file_line_number = 0;
file_indx_number = 1; file_indx_number = 0;
% Simulations. % Simulations.
while iteration <= NumberOfSimulations while iteration < NumberOfSimulations
loop_indx = loop_indx+1; loop_indx = loop_indx+1;
if ( (file_line_number-1)==TableOfInformations(file_indx_number,2))
save([ PriorDirectoryName '/prior_draws' int2str(file_indx_number) '.mat' ],'pdraws');
file_line_number = 1;
file_indx_number = file_indx_number + 1;
pdraws = cell(TableOfInformations(file_indx_number,2),drsave+1) ;
end
params = prior_draw(); params = prior_draw();
set_all_parameters(params); set_all_parameters(params);
[dr,INFO] = resol(oo_.steady_state,work); [dr,INFO] = resol(oo_.steady_state,work);
switch INFO(1) switch INFO(1)
case 0 case 0
file_line_number = file_line_number + 1 ;
iteration = iteration + 1;
pdraws(file_line_number,1) = {params}; pdraws(file_line_number,1) = {params};
if drsave if drsave
pdraws(file_line_number,2) = {dr}; pdraws(file_line_number,2) = {dr};
end end
iteration = iteration+1;
file_line_number = file_line_number+1;
[sampled_prior_expectation,sampled_prior_covariance] = ... [sampled_prior_expectation,sampled_prior_covariance] = ...
recursive_prior_moments(sampled_prior_expectation,sampled_prior_covariance,params,iteration) ; recursive_prior_moments(sampled_prior_expectation,sampled_prior_covariance,params,iteration);
case 1 case 1
count_static_undefined = count_static_undefined + 1; count_static_undefined = count_static_undefined + 1;
case 2 case 2
@ -118,6 +109,14 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_)
otherwise otherwise
count_unknown_problem = count_unknown_problem + 1 ; count_unknown_problem = count_unknown_problem + 1 ;
end end
if ( file_line_number==TableOfInformations(file_indx_number+1,2) )
file_indx_number = file_indx_number + 1;
save([ PriorDirectoryName '/prior_draws' int2str(file_indx_number) '.mat' ],'pdraws');
if file_indx_number<NumberOfFiles
pdraws = cell(TableOfInformations(file_indx_number+1,2),drsave+1);
end
file_line_number = 0;
end
end end
% Get informations about BK conditions and other things... % Get informations about BK conditions and other things...
@ -141,8 +140,7 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_)
results.prior.mean = sampled_prior_expectation; results.prior.mean = sampled_prior_expectation;
results.prior.variance = sampled_prior_covariance; results.prior.variance = sampled_prior_covariance;
results.prior.mass = 1-results.garbage_share; results.prior.mass = 1-results.garbage_share;
function [mu,sigma] = recursive_prior_moments(m0,s0,newobs,iter) function [mu,sigma] = recursive_prior_moments(m0,s0,newobs,iter)
% Recursive estimation of order one and two moments (expectation and % Recursive estimation of order one and two moments (expectation and
% covariance matrix). newobs should be a row vector. I do not use the % covariance matrix). newobs should be a row vector. I do not use the