From fb3899018e2f6342caf9b4a395f542da76f7136a Mon Sep 17 00:00:00 2001 From: adjemian Date: Mon, 1 Sep 2008 11:44:26 +0000 Subject: [PATCH] + Added trace plots. + Correction of a bug related to ferhat's recent commits (options_.model_mode has to be declared somewhere). git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2013 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/GetAllPosteriorDraws.m | 75 ++++++++++++++------- matlab/global_initialization.m | 3 + matlab/name2index.m | 117 +++++++++++++++++++++++++++++++++ matlab/trace_plot.m | 109 ++++++++++++++++++++++++++++++ 4 files changed, 280 insertions(+), 24 deletions(-) create mode 100644 matlab/name2index.m create mode 100644 matlab/trace_plot.m diff --git a/matlab/GetAllPosteriorDraws.m b/matlab/GetAllPosteriorDraws.m index 4d0f27ead..0b2f69e4b 100644 --- a/matlab/GetAllPosteriorDraws.m +++ b/matlab/GetAllPosteriorDraws.m @@ -1,4 +1,4 @@ -function Draws = GetAllPosteriorDraws(column,FirstMhFile,FirstLine,TotalNumberOfMhFile,NumberOfDraws) +function Draws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, blck) % function Draws = GetAllPosteriorDraws(column,FirstMhFile,FirstLine,TotalNumberOfMhFile,NumberOfDraws) % Gets all posterior draws @@ -35,30 +35,57 @@ function Draws = GetAllPosteriorDraws(column,FirstMhFile,FirstLine,TotalNumberOf global M_ options_ -nblck = options_.mh_nblck; -iline = FirstLine; +nblck = options_.mh_nblck; + +iline = FirstLine; linee = 1; DirectoryName = CheckPath('metropolis'); -Draws = zeros(NumberOfDraws*nblck,1); -logpo = zeros(NumberOfDraws*nblck,1); -ipost=0; -if column<=0, - column=1; - ipost=1; -end -iline0=iline; -for blck = 1:nblck - iline=iline0; - for file = FirstMhFile:TotalNumberOfMhFile - load([DirectoryName '/' M_.fname '_mh' int2str(file) '_blck' int2str(blck)],'x2','logpo2') - NumberOfLines = size(x2(iline:end,:),1); - Draws(linee:linee+NumberOfLines-1) = x2(iline:end,column); - logpo(linee:linee+NumberOfLines-1) = logpo2(iline:end); - linee = linee+NumberOfLines; - iline = 1; - end -end -if ipost, - Draws=logpo; +if nblck>1 && nargin<6 + Draws = zeros(NumberOfDraws*nblck,1); + iline0=iline; + if column>0 + for blck = 1:nblck + iline=iline0; + for file = FirstMhFile:TotalNumberOfMhFile + load([DirectoryName '/' M_.fname '_mh' int2str(file) '_blck' int2str(blck)],'x2') + NumberOfLines = size(x2(iline:end,:),1); + Draws(linee:linee+NumberOfLines-1) = x2(iline:end,column); + linee = linee+NumberOfLines; + iline = 1; + end + end + else + for blck = 1:nblck + iline=iline0; + for file = FirstMhFile:TotalNumberOfMhFile + load([DirectoryName '/' M_.fname '_mh' int2str(file) '_blck' int2str(blck)],'logpo2') + NumberOfLines = size(logpo2(iline:end),1); + Draws(linee:linee+NumberOfLines-1) = logpo2(iline:end); + linee = linee+NumberOfLines; + iline = 1; + end + end + end +else + if nblck==1 + blck=1; + end + if column>0 + for file = FirstMhFile:TotalNumberOfMhFile + load([DirectoryName '/' M_.fname '_mh' int2str(file) '_blck' int2str(blck)],'x2') + NumberOfLines = size(x2(iline:end,:),1); + Draws(linee:linee+NumberOfLines-1) = x2(iline:end,column); + linee = linee+NumberOfLines; + iline = 1; + end + else + for file = FirstMhFile:TotalNumberOfMhFile + load([DirectoryName '/' M_.fname '_mh' int2str(file) '_blck' int2str(blck)],'logpo2') + NumberOfLines = size(logpo2(iline:end,:),1); + Draws(linee:linee+NumberOfLines-1) = logpo2(iline:end); + linee = linee+NumberOfLines; + iline = 1; + end + end end \ No newline at end of file diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 8d77460fb..39fd687c3 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -114,6 +114,7 @@ function global_initialization() options_.replic = 50; options_.drop = 100; options_.simul_algo = 0; + options_.model_mode = 0; % if mjdgges.dll (or .mexw32 or ....) doesn't exit, matlab/qz is added to the path. % There exists now qz/mjdgges.m that contains the calls to the old Sims code % Hence, if mjdgges.m is visible exist(...)==2, @@ -180,6 +181,8 @@ function global_initialization() options_.posterior_sampling_method = 'random_walk_metropolis_hastings'; options_.proposal_distribution = 'rand_multivariate_normal'; options_.student_degrees_of_freedom = 3; + options_.trace_plot_ma = 50; + % Misc options_.conf_sig = 0.6; diff --git a/matlab/name2index.m b/matlab/name2index.m new file mode 100644 index 000000000..2eeda95d2 --- /dev/null +++ b/matlab/name2index.m @@ -0,0 +1,117 @@ +function i = name2index(options_, M_, estim_params_, type, name1, name2 ) +% Returns the index associated to an estimated object (deep parameter, +% variance of a structural shock or measurement error, covariance between +% two structural shocks, covariance between two measurement errors). +% +% INPUTS: +% options_ [structure] Dynare structure. +% M_ [structure] Dynare structure (related to model definition). +% estim_params_ [structure] Dynare structure (related to estimation). +% type [string] 'DeepParameter', 'MeasurementError' (for measurement equation error) or 'StructuralShock' (for structural shock). +% name1 [string] +% name2 [string] +% OUTPUTS +% i [integer] column index (in x2, an array of mh draws) associated to name1[,name2]. +% +% SPECIAL REQUIREMENTS +% none + +% Copyright (C) 2008 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 . + + nvx = estim_params_.nvx; + nvn = estim_params_.nvn; + ncx = estim_params_.ncx; + ncn = estim_params_.ncn; + npa = estim_params_.np ; + nnn = nvx+nvn+ncx+ncn+npa; + + i = []; + + if strcmpi(type,'DeepParameter') + i = nvx + nvn + ncx + ncn + ... + strmatch(name1,M_.param_names(estim_params_.param_vals(:,1),:),'exact'); + if nargin>5 + disp('The last input argument is useless!') + end + if isempty(i) + disp([name1 ' is not an estimated deep parameter!']) + end + return + end + + if strcmpi(type,'StructuralShock') + if nargin<6% Covariance matrix diagonal term. + i = strmatch(name1,M_.exo_names(estim_params_.var_exo(:,1),:),'exact'); + if isempty(i) + disp(['The standard deviation of ' name1 ' is not an estimated parameter!']) + return + end + else% Covariance matrix off-diagonal term + offset = nvx+nvn; + try + list_of_structural_shocks = [ M_.exo_names(estim_params_.corrx(:,1),:) , M_.exo_names(estim_params_.corrx(:,2),:) ]; + k1 = strmatch(name1,list_of_structural_shocks(:,1),'exact'); + k2 = strmatch(name2,list_of_structural_shocks(:,2),'exact'); + i = offset+intersect(k1,k2); + if isempty(i) + k1 = strmatch(name1,list_of_structural_shocks(:,2),'exact'); + k2 = strmatch(name2,list_of_structural_shocks(:,1),'exact'); + i = offset+intersect(k1,k2); + end + if isempty(i) + if isempty(i) + disp(['The correlation between' name1 ' and ' name2 ' is not an estimated parameter!']) + return + end + end + catch + disp(['Off diagonal terms of the covariance matrix are not estimated (state equation)']) + end + end + end + + if strcmpi(type,'MeasurementError') + if nargin<6% Covariance matrix diagonal term + i = nvx + strmatch(name1,M_.endo_names(estim_params_.var_endo(:,1),:),'exact'); + if isempty(i) + disp(['The standard deviation of the measurement error on ' name1 ' is not an estimated parameter!']) + return + end + else% Covariance matrix off-diagonal term + offset = nvx+nvn+ncx; + try + list_of_measurement_errors = [ M_.endo_names(estim_params_.corrn(:,1),:) , M_.endo_names(estim_params_.corrn(:,2),:) ]; + k1 = strmatch(name1,list_of_measurement_errors(:,1),'exact'); + k2 = strmatch(name2,list_of_measurement_errors(:,2),'exact'); + i = offset+intersect(k1,k2); + if isempty(i) + k1 = strmatch(name1,list_of_measurement_errors(:,2),'exact'); + k2 = strmatch(name2,list_of_measurement_errors(:,1),'exact'); + i = offset+intersect(k1,k2); + end + if isempty(i) + if isempty(i) + disp(['The correlation between the measurement errors on ' name1 ' and ' name2 ' is not an estimated parameter!']) + return + end + end + catch + disp('Off diagonal terms of the covariance matrix are not estimated (measurement equation)') + end + end + end \ No newline at end of file diff --git a/matlab/trace_plot.m b/matlab/trace_plot.m new file mode 100644 index 000000000..78216c5c8 --- /dev/null +++ b/matlab/trace_plot.m @@ -0,0 +1,109 @@ +function trace_plot(options_,M_,estim_params_,type,blck,name1,name2) +% This function build trace plot for the metropolis hastings draws. +% +% +% INPUTS +% +% name +% +% OUTPUTS +% density [double] density +% +% SPECIAL REQUIREMENTS + +% Copyright (C) 2003-2008 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 . + +% Cet the column index: +if nargin<7 + column = name2index(options_, M_, estim_params_, type, name1); +else + column = name2index(options_, M_, estim_params_, type, name1, name2); +end + +if isempty(column) + return +end + +% Get informations about the posterior draws: +DirectoryName = CheckPath('metropolis'); +try + load([DirectoryName '/' M_.fname '_mh_history.mat']); +catch + disp(['trace_plot:: I can''t find ' M_.fname '_results.mat !']) + disp(['trace_plot:: Did you run a metropolis?']) + return +end + +FirstMhFile = 1; +FirstLine = 1; +TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); +TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); +clear record; + +% Get all the posterior draws: +PosteriorDraws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, blck); + + +% Plot the posterior draws: + +if strcmpi(type,'DeepParameter') + TYPE = 'parameter '; +elseif strcmpi(type,'StructuralShock') + if nargin<7 + TYPE = 'the standard deviation of structural shock '; + else + TYPE = 'the correlation between structural shocks '; + end +elseif strcmpi(type,'MeasurementError') + if nargin<7 + TYPE = 'the standard deviation of measurement error '; + else + TYPE = 'the correlation between measurement errors '; + end +end + +if nargin<7 + FigureName = ['trace plot for ' TYPE name1]; +else + FigureName = ['trace plot for ' TYPE name1 ' and ' name2]; +end + +if options_.mh_nblck>1 + FigureName = [ FigureName , ' (block number' int2str(blck) ').']; +end + + +figure('Name',FigureName) +plot(1:TotalNumberOfMhDraws,PosteriorDraws,'Color',[.8 .8 .8]); + + +% Compute the moving average of the posterior draws: + +N = options_.trace_plot_ma; +MovingAverage = NaN(TotalNumberOfMhDraws,1); +first = N+1; +last = TotalNumberOfMhDraws-N; + +for t=first:last + MovingAverage(t) = mean(PosteriorDraws(t-N:t+N)); +end + +hold on +plot(1:TotalNumberOfMhDraws,MovingAverage,'-k','linewidth',2) +hold off +axis tight \ No newline at end of file