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