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

Cosmetic changes to various estimation routines

kalman-mex
Johannes Pfeifer 2023-07-13 14:05:52 -04:00
parent e6c43c2a29
commit 8532d6abd7
9 changed files with 65 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
matlab/GetPosteriorMeanVariance.m
View File

@ -1,6 +1,20 @@
function [mean,variance] = GetPosteriorMeanVariance(M,drop) function [mean,variance] = GetPosteriorMeanVariance(M_,drop)
%function [mean,variance] = GetPosteriorMeanVariance(M,drop)
% Computes the posterior mean and variance
% (+updates of oo_ & TeX output).
%
% INPUTS
% M_ [structure] Dynare model structure
% drop [double] share of draws to drop
%
% OUTPUTS
% mean [double] mean parameter vector
% variance [double] variance
%
% SPECIAL REQUIREMENTS
% None.
% Copyright © 2012-2017 Dynare Team % Copyright © 2012-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -17,8 +31,8 @@ function [mean,variance] = GetPosteriorMeanVariance(M,drop)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
MetropolisFolder = CheckPath('metropolis',M.dname); MetropolisFolder = CheckPath('metropolis',M_.dname);
FileName = M.fname; FileName = M_.fname;
BaseName = [MetropolisFolder filesep FileName]; BaseName = [MetropolisFolder filesep FileName];
record=load_last_mh_history_file(MetropolisFolder, FileName); record=load_last_mh_history_file(MetropolisFolder, FileName);
NbrDraws = sum(record.MhDraws(:,1)); NbrDraws = sum(record.MhDraws(:,1));

24
matlab/GetPosteriorParametersStatistics.m
View File

@ -1,17 +1,18 @@
function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames) function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames)
% function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames)
% This function prints and saves posterior estimates after the mcmc % This function prints and saves posterior estimates after the mcmc
% (+updates of oo_ & TeX output). % (+updates of oo_ & TeX output).
% %
% INPUTS % INPUTS
% estim_params_ [structure] % estim_params_ [structure] Dynare estimation parameter structure
% M_ [structure] % M_ [structure] Dynare model structure
% options_ [structure] % options_ [structure] Dynare options structure
% bayestopt_ [structure] % bayestopt_ [structure] Dynare structure describing priors
% oo_ [structure] % oo_ [structure] Dynare results structure
% pnames [cell] cell of strings, names of the prior shapes available % pnames [cell] cell of strings, names of the prior shapes available
% %
% OUTPUTS % OUTPUTS
% oo_ [structure] % oo_ [structure] Dynare results structure
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% None. % None.
@ -50,7 +51,6 @@ FileName = M_.fname;
record=load_last_mh_history_file(MetropolisFolder,FileName); record=load_last_mh_history_file(MetropolisFolder,FileName);
FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
@ -68,12 +68,10 @@ skipline(2)
disp('ESTIMATION RESULTS') disp('ESTIMATION RESULTS')
skipline() skipline()
try if ~isfield(oo_,'MarginalDensity') || ~isfield(oo_.MarginalDensity,'ModifiedHarmonicMean')
disp(sprintf('Log data density is %f.', oo_.MarginalDensity.ModifiedHarmonicMean)) [~,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
catch
[marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
disp(sprintf('Log data density is %f.', oo_.MarginalDensity.ModifiedHarmonicMean))
end end
fprintf('\nLog data density is %f.\n', oo_.MarginalDensity.ModifiedHarmonicMean);
num_draws=NumberOfDraws*mh_nblck; num_draws=NumberOfDraws*mh_nblck;
hpd_draws = round((1-options_.mh_conf_sig)*num_draws); hpd_draws = round((1-options_.mh_conf_sig)*num_draws);

18
matlab/convergence_diagnostics/McMCDiagnostics.m
View File

@ -1,11 +1,12 @@
function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_) function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_)
% function McMCDiagnostics % function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_)
% Computes convergence tests % Computes convergence tests
% %
% INPUTS % INPUTS
% options_ [structure] % options_ [structure] Dynare options structure
% estim_params_ [structure] % estim_params_ [structure] Dynare estimation parameter structure
% M_ [structure] % M_ [structure] Dynare model structure
% oo_ [structure] Dynare results structure
% %
% OUTPUTS % OUTPUTS
% oo_ [structure] % oo_ [structure]
@ -58,7 +59,6 @@ if issue_an_error_message
end end
% compute inefficiency factor % compute inefficiency factor
FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
@ -68,6 +68,7 @@ NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws
param_name = {}; param_name = {};
param_name_tex = {}; param_name_tex = {};
Ifac=NaN(nblck,npar);
for jj = 1:npar for jj = 1:npar
if options_.TeX if options_.TeX
[par_name_temp, par_name_tex_temp] = get_the_name(jj, options_.TeX, M_,estim_params_, options_); [par_name_temp, par_name_tex_temp] = get_the_name(jj, options_.TeX, M_,estim_params_, options_);
@ -107,8 +108,6 @@ update_last_mh_history_file(MetropolisFolder, ModelName, record);
PastDraws = sum(record.MhDraws,1); PastDraws = sum(record.MhDraws,1);
LastFileNumber = PastDraws(2);
LastLineNumber = record.MhDraws(end,3);
NumberOfDraws = PastDraws(1); NumberOfDraws = PastDraws(1);
if NumberOfDraws<=2000 if NumberOfDraws<=2000
@ -217,8 +216,6 @@ ALPHA = 0.2; % increase too much with the number
time = 1:NumberOfDraws; time = 1:NumberOfDraws;
xx = Origin:StepSize:NumberOfDraws; xx = Origin:StepSize:NumberOfDraws;
NumberOfLines = length(xx); NumberOfLines = length(xx);
tmp = zeros(NumberOfDraws*nblck,3);
UDIAG = zeros(NumberOfLines,6,npar);
if NumberOfDraws < Origin if NumberOfDraws < Origin
disp('Estimation::mcmc::diagnostics: The number of simulations is too small to compute the MCMC convergence diagnostics.') disp('Estimation::mcmc::diagnostics: The number of simulations is too small to compute the MCMC convergence diagnostics.')
@ -256,7 +253,6 @@ if isnumeric(options_.parallel)
clear fout clear fout
% Parallel execution! % Parallel execution!
else else
ModelName = ModelName;
if ~isempty(M_.bvar) if ~isempty(M_.bvar)
ModelName = [ModelName '_bvar']; ModelName = [ModelName '_bvar'];
end end
@ -391,7 +387,7 @@ if reste
fprintf(fidTeX,'\\label{Fig:UnivariateDiagnostics:%s}\n',int2str(pages+1)); fprintf(fidTeX,'\\label{Fig:UnivariateDiagnostics:%s}\n',int2str(pages+1));
fprintf(fidTeX,'\\end{figure}\n'); fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,'\n'); fprintf(fidTeX,'\n');
fprintf(fidTeX,'% End Of TeX file.'); fprintf(fidTeX,'%% End Of TeX file.');
fclose(fidTeX); fclose(fidTeX);
end end
end % if reste > 0 end % if reste > 0

2
matlab/get_name_of_the_last_mh_file.m
View File

@ -1,4 +1,5 @@
function [mhname,info] = get_name_of_the_last_mh_file(M_) function [mhname,info] = get_name_of_the_last_mh_file(M_)
% function [mhname,info] = get_name_of_the_last_mh_file(M_)
% This function returns the name of the last mh file and test if the metropolis was completed. % This function returns the name of the last mh file and test if the metropolis was completed.
% %
% INPUTS % INPUTS
@ -28,7 +29,6 @@ function [mhname,info] = get_name_of_the_last_mh_file(M_)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
mhname = [];
info = 1; info = 1;
MetropolisFolder = CheckPath('metropolis',M_.dname); MetropolisFolder = CheckPath('metropolis',M_.dname);

17
matlab/marginal_density.m
View File

@ -1,16 +1,16 @@
function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_) function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_)
% function marginal = marginal_density() % function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_)
% Computes the marginal density % Computes the marginal density
% %
% INPUTS % INPUTS
% options_ [structure] % options_ [structure] Dynare options structure
% estim_params_ [structure] % estim_params_ [structure] Dynare estimation parameter structure
% M_ [structure] % M_ [structure] Dynare model structure
% oo_ [structure] % oo_ [structure] Dynare results structure
% %
% OUTPUTS % OUTPUTS
% marginal: [double] marginal density (modified harmonic mean) % marginal: [double] marginal density (modified harmonic mean)
% oo_ [structure] % oo_ [structure] Dynare results structure
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
@ -42,10 +42,9 @@ record=load_last_mh_history_file(MetropolisFolder, ModelName);
[nblck, npar] = size(record.LastParameters); [nblck, npar] = size(record.LastParameters);
FirstMhFile = record.KeepedDraws.FirstMhFile; FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine; FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
TODROP = floor(options_.mh_drop*TotalNumberOfMhDraws); TODROP = floor(options_.mh_drop*TotalNumberOfMhDraws);
fprintf('Estimation::marginal density: I''m computing the posterior mean and covariance... '); fprintf('Estimation::marginal density: I''m computing the posterior mean and covariance... ');

4
matlab/metropolis_draw.m
View File

@ -1,5 +1,5 @@
function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_) function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_)
% function [xparams, logpost]=metropolis_draw(init) % function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_)
% Builds draws from metropolis % Builds draws from metropolis
% %
% INPUTS: % INPUTS:
@ -62,8 +62,6 @@ if init
record=load_last_mh_history_file(MetropolisFolder, FileName); record=load_last_mh_history_file(MetropolisFolder, FileName);
FirstMhFile = record.KeepedDraws.FirstMhFile; FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
LastMhFile = TotalNumberOfMhFiles;
TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); %number of parameters plus posterior plus ? MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); %number of parameters plus posterior plus ?

29
matlab/posterior_sampler.m
View File

@ -3,21 +3,20 @@ function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_boun
% Random Walk Metropolis-Hastings algorithm. % Random Walk Metropolis-Hastings algorithm.
% %
% INPUTS % INPUTS
% o TargetFun [char] string specifying the name of the objective % o TargetFun [char] string specifying the name of the objective
% function (posterior kernel). % function (posterior kernel).
% o ProposalFun [char] string specifying the name of the proposal % o ProposalFun [char] string specifying the name of the proposal
% density % density
% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values). % o xparam1 [double] (p*1) vector of parameters to be estimated (initial values).
% o sampler_options structure % o sampler_options structure
% .invhess [double] (p*p) matrix, posterior covariance matrix (at the mode). % o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters.
% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters. % o dataset_ [structure] data structure
% o dataset_ data structure % o dataset_info [structure] dataset info structure
% o dataset_info dataset info structure % o options_ [structure] options structure
% o options_ options structure % o M_ [structure] model structure
% o M_ model structure % o estim_params_ [structure] estimated parameters structure
% o estim_params_ estimated parameters structure % o bayestopt_ [structure] prior specification structure
% o bayestopt_ estimation options structure % o oo_ [structure] output structure
% o oo_ outputs structure
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% None. % None.

9
matlab/posterior_sampler_initialization.m
View File

@ -58,16 +58,9 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npa
%Initialize outputs %Initialize outputs
ix2 = []; ix2 = [];
ilogpo2 = []; ilogpo2 = [];
ModelName = [];
MetropolisFolder = [];
FirstBlock = []; FirstBlock = [];
FirstLine = []; FirstLine = [];
npar = [];
NumberOfBlocks = [];
nruns = [];
NewFile = []; NewFile = [];
MAX_nruns = [];
d = [];
ModelName = M_.fname; ModelName = M_.fname;
if ~isempty(M_.bvar) if ~isempty(M_.bvar)
@ -171,7 +164,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
else else
bayestopt0 = load([PreviousFolder0 filesep 'prior' filesep 'definition.mat']); bayestopt0 = load([PreviousFolder0 filesep 'prior' filesep 'definition.mat']);
end end
[common_parameters,IA,IB] = intersect(bayestopt_.name,bayestopt0.bayestopt_.name); [~,IA,IB] = intersect(bayestopt_.name,bayestopt0.bayestopt_.name);
new_estimated_parameters = ~ismember(bayestopt_.name,bayestopt0.bayestopt_.name); new_estimated_parameters = ~ismember(bayestopt_.name,bayestopt0.bayestopt_.name);
ix2 = zeros(NumberOfBlocks,npar); ix2 = zeros(NumberOfBlocks,npar);
ilogpo2 = zeros(NumberOfBlocks,1); ilogpo2 = zeros(NumberOfBlocks,1);

10
matlab/prior_posterior_statistics.m
View File

@ -1,6 +1,5 @@
function prior_posterior_statistics(type,dataset,dataset_info) function prior_posterior_statistics(type,dataset,dataset_info)
% function prior_posterior_statistics(type,dataset,dataset_info)
% function prior_posterior_statistics(type,dataset)
% Computes Monte Carlo filter smoother and forecasts % Computes Monte Carlo filter smoother and forecasts
% %
% INPUTS % INPUTS
@ -8,6 +7,7 @@ function prior_posterior_statistics(type,dataset,dataset_info)
% prior % prior
% gsa % gsa
% dataset: data structure % dataset: data structure
% dataset_info: dataset structure
% %
% OUTPUTS % OUTPUTS
% none % none
@ -208,13 +208,11 @@ localVars.ifil=ifil;
localVars.DirectoryName = DirectoryName; localVars.DirectoryName = DirectoryName;
if strcmpi(type,'posterior') if strcmpi(type,'posterior')
BaseName = [DirectoryName filesep M_.fname];
record=load_last_mh_history_file(DirectoryName, M_.fname); record=load_last_mh_history_file(DirectoryName, M_.fname);
[nblck, npar] = size(record.LastParameters); [nblck, npar] = size(record.LastParameters);
FirstMhFile = record.KeepedDraws.FirstMhFile; FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
LastMhFile = TotalNumberOfMhFiles;
TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
mh_nblck = options_.mh_nblck; mh_nblck = options_.mh_nblck;
@ -318,9 +316,9 @@ if strcmpi(type,'gsa')
return return
end end
if ~isnumeric(options_.parallel), if ~isnumeric(options_.parallel)
leaveSlaveOpen = options_.parallel_info.leaveSlaveOpen; leaveSlaveOpen = options_.parallel_info.leaveSlaveOpen;
if options_.parallel_info.leaveSlaveOpen == 0, if options_.parallel_info.leaveSlaveOpen == 0
% Commenting for testing!!! % Commenting for testing!!!
options_.parallel_info.leaveSlaveOpen = 1; % Force locally to leave open remote matlab sessions (repeated pm3 calls) options_.parallel_info.leaveSlaveOpen = 1; % Force locally to leave open remote matlab sessions (repeated pm3 calls)
end end