From cd5e888f08d4c912dbc05c91a3cc0ddb2a007594 Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Tue, 19 Mar 2013 21:12:49 +0100
Subject: [PATCH] Fix several bugs related to estimated measurement errors
1. The first call to set_prior overwrote the first column of
estim_params_.var_endo storing the position of the variable with
measurement error in M_.endo_names with the position in
options_.var_obs. All subsequent calls to set_prior then lead to
crashes.
2. At the same time, for correlations of ME, the first column of
estim_params_.corrn still stored the position of the variable with
measurement error in M_.endo_names. But subsequent calls to it were done
as if it stored the position in options_.var_obs
I introduced two new variables in estim_params_ storing the respective
positions in var_obs so as to not necessitate changes in the
preprocessors.
3. For cases of calibrated measurement error correlations, the
covariance matrix was not updated.
4. Fixing a lot of smaller bugs related to measurement errors, including
some copy and paste errors
-
(cherry picked from commit 7518072e7762c54e39d75c2f6e1e6261f416efa4)
---
matlab/GetPosteriorParametersStatistics.m | 10 ++++----
matlab/PlotPosteriorDistributions.m | 4 +--
matlab/dynare_estimation_1.m | 31 ++++++++++++++++++-----
matlab/get_posterior_parameters.m | 8 +++---
matlab/get_the_name.m | 8 +++---
matlab/global_initialization.m | 1 +
matlab/prior_posterior_statistics.m | 10 ++++----
matlab/row_header_width.m | 4 +--
matlab/set_all_parameters.m | 16 +++++++-----
matlab/set_prior.m | 29 +++++++++++++--------
10 files changed, 77 insertions(+), 44 deletions(-)
diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m
index ef706b584..cb6ba4219 100644
--- a/matlab/GetPosteriorParametersStatistics.m
+++ b/matlab/GetPosteriorParametersStatistics.m
@@ -15,7 +15,7 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bay
% SPECIAL REQUIREMENTS
% None.
-% Copyright (C) 2006-2012 Dynare Team
+% Copyright (C) 2006-2013 Dynare Team
%
% This file is part of Dynare.
%
@@ -163,7 +163,7 @@ end
if nvn
type = 'measurement_errors_std';
if TeX
- fid = TeXBegin(OutputDirectoryName,M_.fname,3,'standard deviation of measurement errors')
+ fid = TeXBegin(OutputDirectoryName,M_.fname,3,'standard deviation of measurement errors');
end
disp(' ')
disp('standard deviation of measurement errors')
@@ -174,17 +174,17 @@ if nvn
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
posterior_moments(Draws,1,options_.mh_conf_sig);
- name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
+ name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:));
oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
else
try
- name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
+ name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:));
[post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
catch
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
posterior_moments(Draws,1,options_.mh_conf_sig);
- name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
+ name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:));
oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
end
end
diff --git a/matlab/PlotPosteriorDistributions.m b/matlab/PlotPosteriorDistributions.m
index ba6510eea..5a622272b 100644
--- a/matlab/PlotPosteriorDistributions.m
+++ b/matlab/PlotPosteriorDistributions.m
@@ -16,7 +16,7 @@ function oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt
% SPECIAL REQUIREMENTS
% none
-% Copyright (C) 2005-2012 Dynare Team
+% Copyright (C) 2005-2013 Dynare Team
%
% This file is part of Dynare.
%
@@ -89,7 +89,7 @@ for i=1:npar
eval(['pmod = oo_.posterior_mode.shocks_std.' name ';'])
end
elseif i <= nvx+nvn
- name = deblank(options_.varobs(estim_params_.var_endo(i-nvx,1),:));
+ name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i-nvx,1),:));
eval(['x1 = oo_.posterior_density.measurement_errors_std.' name '(:,1);'])
eval(['f1 = oo_.posterior_density.measurement_errors_std.' name '(:,2);'])
eval(['oo_.prior_density.mearsurement_errors_std.' name '(:,1) = x2;'])
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 1b4338b78..d5ea73a35 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -77,6 +77,24 @@ if ~isequal(estim_params_.ncx,nnz(tril(M_.Sigma_e,-1)))
end
end
+M_.H_is_diagonal = 1;
+if estim_params_.ncn || ~isequal(nnz(M_.H),length(M_.H))
+ M_.H_is_diagonal = 0;
+end
+
+% Set the correlation matrix of measurement errors if necessary.
+if ~isequal(estim_params_.ncn,nnz(tril(M_.H,-1)))
+ M_.Correlation_matrix_ME = diag(1./sqrt(diag(M_.H)))*M_.H*diag(1./sqrt(diag(M_.H)));
+ % Remove NaNs appearing because of variances calibrated to zero.
+ if any(isnan(M_.Correlation_matrix_ME))
+ zero_variance_idx = find(~diag(M_.H));
+ for i=1:length(zero_variance_idx)
+ M_.Correlation_matrix_ME(zero_variance_idx(i),:) = 0;
+ M_.Correlation_matrix_ME(:,zero_variance_idx(i)) = 0;
+ end
+ end
+end
+
data = dataset_.data;
rawdata = dataset_.rawdata;
data_index = dataset_.missing.aindex;
@@ -543,7 +561,7 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
disp(tit1)
ip = nvx+1;
for i=1:nvn
- name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
+ name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:));
disp(sprintf('%-*s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
header_width,name,bayestopt_.p1(ip), ...
xparam1(ip),stdh(ip),tstath(ip), ...
@@ -642,7 +660,7 @@ elseif ~any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
disp(tit1)
ip = nvx+1;
for i=1:nvn
- name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
+ name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:));
disp(sprintf('%-*s %8.4f %7.4f %7.4f',header_width,name,xparam1(ip),stdh(ip),tstath(ip)))
eval(['oo_.mle_mode.measurement_errors_std.' name ' = xparam1(ip);']);
eval(['oo_.mle_std.measurement_errors_std.' name ' = stdh(ip);']);
@@ -788,7 +806,7 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
fprintf(fidTeX,'\\hline \\hline \\endlastfoot \n');
ip = nvx+1;
for i=1:nvn
- idx = strmatch(options_.varobs(estim_params_.var_endo(i,1),:),M_.endo_names);
+ idx = strmatch(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:),M_.endo_names);
fprintf(fidTeX,'$%s$ & %4s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',...
deblank(M_.endo_names_tex(idx,:)), ...
deblank(pnames(bayestopt_.pshape(ip)+1,:)), ...
@@ -1089,7 +1107,8 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
hh = dyn_figure(options_,'Name','Smoothed observation errors');
NAMES = [];
if options_.TeX, TeXNAMES = []; end
- for i=1:min(nstar,number_of_plots_to_draw-(nbplt-1)*nstar)
+ nstar0=min(nstar,number_of_plots_to_draw-(nbplt-1)*nstar);
+ for i=1:nstar0
k = (plt-1)*nstar+i;
subplot(nr,nc,i);
plot([1 gend],[0 0],'-r','linewidth',.5)
@@ -1120,7 +1139,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
dyn_saveas(hh,[M_.fname '_SmoothedObservationErrors' int2str(plt)],options_);
if options_.TeX
fprintf(fidTeX,'\\begin{figure}[H]\n');
- for jj = 1:nstar
+ for jj = 1:nstar0
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
end
fprintf(fidTeX,'\\centering \n');
@@ -1157,7 +1176,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
for i=1:nstar0,
k = (plt-1)*nstar+i;
subplot(nr,nc,i);
- plot(1:gend,yf(k,:),'--r','linewidth',1)
+ plot(1:gend,yf(k,:),'-r','linewidth',1)
hold on
plot(1:gend,rawdata(:,k),'--k','linewidth',1)
hold off
diff --git a/matlab/get_posterior_parameters.m b/matlab/get_posterior_parameters.m
index dc987f610..8c4a9c8eb 100644
--- a/matlab/get_posterior_parameters.m
+++ b/matlab/get_posterior_parameters.m
@@ -12,7 +12,7 @@ function xparam = get_posterior_parameters(type)
% SPECIAL REQUIREMENTS
% None.
-% Copyright (C) 2006-2009 Dynare Team
+% Copyright (C) 2006-2013 Dynare Team
%
% This file is part of Dynare.
%
@@ -49,7 +49,7 @@ for i=1:nvx
end
for i=1:nvn
- k1 = estim_params_.var_endo(i,1);
+ k1 = estim_params_.nvn_observable_correspondence(i,1);
name1 = deblank(options_.varobs(k1,:));
xparam(m) = eval(['oo_.posterior_' type '.measurement_errors_std.' name1]);
m = m+1;
@@ -67,8 +67,8 @@ for i=1:ncx
end
for i=1:ncn
- k1 = estim_params_.corrn(i,1);
- k2 = estim_params_.corrn(i,2);
+ k1 = estim_params_.corrn_observable_correspondence(i,1);
+ k2 = estim_params_.corrn_observable_correspondence(i,2);
name1 = deblank(options_.varobs(k1,:));
name2 = deblank(options_.varobs(k2,:));
xparam(m) = eval(['oo_.posterior_' type '.measurement_errors_corr.' name1 '_' name2]);
diff --git a/matlab/get_the_name.m b/matlab/get_the_name.m
index 9455b7181..bd6172a22 100644
--- a/matlab/get_the_name.m
+++ b/matlab/get_the_name.m
@@ -40,7 +40,7 @@ function [nam,texnam] = get_the_name(k,TeX,M_,estim_params_,options_)
%! @end deftypefn
%@eod:
-% Copyright (C) 2004-2011 Dynare Team
+% Copyright (C) 2004-2013 Dynare Team
%
% This file is part of Dynare.
%
@@ -73,11 +73,11 @@ if k <= nvx
texnam = ['$ SE_{' tname '} $'];
end
elseif k <= (nvx+nvn)
- vname = deblank(options_.varobs(estim_params_.var_endo(k-estim_params_.nvx,1),:));
+ vname = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(k-estim_params_.nvx,1),:));
nam=['SE_EOBS_',vname];
if TeX
tname = deblank(M_.endo_names_tex(estim_params_.var_endo(k-estim_params_.nvx,1),:));
- texnam = ['$ SE_{' tname '} $'];
+ texnam = ['$ EOBS SE_{' tname '} $'];
end
elseif k <= (nvx+nvn+ncx)
jj = k - (nvx+nvn);
@@ -97,7 +97,7 @@ elseif k <= (nvx+nvn+ncx+ncn)
nam=['CC_EOBS_' vname];
if TeX
tname = [deblank(M_.endo_names_tex(k1,:)) ',' deblank(M_.endo_names_tex(k2,:))];
- texnam =['$ CC_{' tname '} $'];
+ texnam =['$ EOBS CC_{' tname '} $'];
end
else
jj = k - (nvx+nvn+ncx+ncn);
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 29783c0e0..ef312e8f8 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -400,6 +400,7 @@ oo_.exo_det_simul = [];
M_.params = [];
M_.endo_histval = [];
M_.Correlation_matrix = [];
+M_.Correlation_matrix_ME = [];
% homotopy
options_.homotopy_mode = 0;
diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m
index 3bb87776b..8eb763508 100644
--- a/matlab/prior_posterior_statistics.m
+++ b/matlab/prior_posterior_statistics.m
@@ -19,7 +19,7 @@ function prior_posterior_statistics(type,dataset)
% See the comments random_walk_metropolis_hastings.m funtion.
-% Copyright (C) 2005-2012 Dynare Team
+% Copyright (C) 2005-2013 Dynare Team
%
% This file is part of Dynare.
%
@@ -298,20 +298,20 @@ end
if options_.filtered_vars
pm3(endo_nbr,gend,ifil(1),B,'Updated Variables',...
- '',varlist,'tit_tex',M_.endo_names,...
+ '',varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'UpdatedVariables',DirectoryName, ...
'_update');
pm3(endo_nbr,gend+1,ifil(4),B,'One step ahead forecast',...
- '',varlist,'tit_tex',M_.endo_names,...
+ '',varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'FilteredVariables',DirectoryName,'_filter_step_ahead');
end
if options_.forecast
pm3(endo_nbr,horizon+maxlag,ifil(6),B,'Forecasted variables (mean)',...
- '',varlist,'tit_tex',M_.endo_names,...
+ '',varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'MeanForecast',DirectoryName,'_forc_mean');
pm3(endo_nbr,horizon+maxlag,ifil(6),B,'Forecasted variables (point)',...
- '',varlist,'tit_tex',M_.endo_names,...
+ '',varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'PointForecast',DirectoryName,'_forc_point');
end
diff --git a/matlab/row_header_width.m b/matlab/row_header_width.m
index 2e07d57f5..92ef63d1e 100644
--- a/matlab/row_header_width.m
+++ b/matlab/row_header_width.m
@@ -13,7 +13,7 @@ function w=row_header_width(M_,estim_params_,bayestopt_)
% SPECIAL REQUIREMENTS
% None.
-% Copyright (C) 2006-2009 Dynare Team
+% Copyright (C) 2006-2013 Dynare Team
%
% This file is part of Dynare.
%
@@ -64,7 +64,7 @@ if ncx
end
end
if ncn
- for i=1:nvn
+ for i=1:ncn
k1 = estim_params_.corrn(i,1);
k2 = estim_params_.corrn(i,2);
w = max(w,length(deblank(M_.endo_names(k1,:)))...
diff --git a/matlab/set_all_parameters.m b/matlab/set_all_parameters.m
index a34fbd9fe..42257ac4e 100644
--- a/matlab/set_all_parameters.m
+++ b/matlab/set_all_parameters.m
@@ -71,9 +71,8 @@ offset = nvx;
% setting measument error variance
if nvn
- var_endo = estim_params.var_endo;
for i=1:nvn
- k = var_endo(i,1);
+ k = estim_params.nvn_observable_correspondence(i,1);
H(k,k) = xparam1(i+offset)^2;
end
end
@@ -100,11 +99,16 @@ end
% update offset
offset = nvx+nvn+ncx;
% setting measurement error covariances
+if ~isempty(M.Correlation_matrix_ME)
+ H = diag(sqrt(diag(H)))*M.Correlation_matrix_ME*diag(sqrt(diag(H)));
+end
if ncn
- corrn = estim_params.corrn;
+ corrn_observable_correspondence = estim_params.corrn_observable_correspondence;
for i=1:ncn
- k1 = corr(i,1);
- k2 = corr(i,2);
+ k1 = corrn_observable_correspondence(i,1);
+ k2 = corrn_observable_correspondence(i,2);
+ M.Correlation_matrix_ME(k1,k2) = xparam1(i+offset);
+ M.Correlation_matrix_ME(k2,k1) = M.Correlation_matrix_ME(k1,k2);
H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
H(k2,k1) = H(k1,k2);
end
@@ -122,6 +126,6 @@ end
if nvx || ncx
M.Sigma_e = Sigma_e;
end
-if nvn
+if nvn || ncn
M.H = H;
end
\ No newline at end of file
diff --git a/matlab/set_prior.m b/matlab/set_prior.m
index 841af7442..f82223ec8 100644
--- a/matlab/set_prior.m
+++ b/matlab/set_prior.m
@@ -18,7 +18,7 @@ function [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params
% SPECIAL REQUIREMENTS
% None
-% Copyright (C) 2003-2011 Dynare Team
+% Copyright (C) 2003-2013 Dynare Team
%
% This file is part of Dynare.
%
@@ -41,11 +41,11 @@ ncx = size(estim_params_.corrx,1);
ncn = size(estim_params_.corrn,1);
np = size(estim_params_.param_vals,1);
-estim_params_.nvx = nvx;
-estim_params_.nvn = nvn;
-estim_params_.ncx = ncx;
-estim_params_.ncn = ncn;
-estim_params_.np = np;
+estim_params_.nvx = nvx; %exogenous shock variances
+estim_params_.nvn = nvn; %endogenous variances, i.e. measurement error
+estim_params_.ncx = ncx; %exogenous shock correlations
+estim_params_.ncn = ncn; % correlation between endogenous variables, i.e. measurement error.
+estim_params_.np = np; % other parameters of the model
xparam1 = [];
ub = [];
@@ -74,6 +74,7 @@ if nvx
bayestopt_.name = cellstr(M_.exo_names(estim_params_.var_exo(:,1),:));
end
if nvn
+ estim_params_.nvn_observable_correspondence=NaN(nvn,1); % stores number of corresponding observable
if isequal(M_.H,0)
nvarobs = size(options_.varobs,1);
M_.H = zeros(nvarobs,nvarobs);
@@ -83,7 +84,7 @@ if nvn
if isempty(obsi_)
error(['The variable ' deblank(M_.endo_names(estim_params_.var_endo(i,1),:)) ' has to be declared as observable since you assume a measurement error on it.'])
end
- estim_params_.var_endo(i,1) = obsi_;
+ estim_params_.nvn_observable_correspondence(i,1)=obsi_;
end
xparam1 = [xparam1; estim_params_.var_endo(:,2)];
ub = [ub; estim_params_.var_endo(:,4)];
@@ -94,7 +95,7 @@ if nvn
bayestopt_.p3 = [ bayestopt_.p3; estim_params_.var_endo(:,8)];
bayestopt_.p4 = [ bayestopt_.p4; estim_params_.var_endo(:,9)];
bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.var_endo(:,10)];
- bayestopt_.name = [ bayestopt_.name; cellstr(options_.varobs(estim_params_.var_endo(:,1),:))];
+ bayestopt_.name = [ bayestopt_.name; cellstr(options_.varobs(estim_params_.nvn_observable_correspondence,:))];
end
if ncx
xparam1 = [xparam1; estim_params_.corrx(:,3)];
@@ -111,6 +112,7 @@ if ncx
repmat(', ',ncx,1) , deblank(M_.exo_names(estim_params_.corrx(:,2),:))])];
end
if ncn
+ estim_params_.corrn_observable_correspondence=NaN(ncn,2);
if isequal(M_.H,0)
nvarobs = size(options_.varobs,1);
M_.H = zeros(nvarobs,nvarobs);
@@ -125,8 +127,15 @@ if ncn
bayestopt_.p4 = [ bayestopt_.p4; estim_params_.corrn(:,10)];
bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.corrn(:,11)];
bayestopt_.name = [bayestopt_.name; cellstr([repmat('corr ',ncn,1) ...
- deblank(M_.exo_names(estim_params_.corrn(:,1),:)) ...
- repmat(', ',ncn,1) , deblank(M_.exo_names(estim_params_.corrn(:,2),:))])];
+ deblank(M_.endo_names(estim_params_.corrn(:,1),:)) ...
+ repmat(', ',ncn,1) , deblank(M_.endo_names(estim_params_.corrn(:,2),:))])];
+ for i=1:ncn
+ k1 = estim_params_.corrn(i,1);
+ k2 = estim_params_.corrn(i,2);
+ obsi1 = strmatch(deblank(M_.endo_names(k1,:)),deblank(options_.varobs),'exact'); %find correspondence to varobs to construct H in set_all_paramters
+ obsi2 = strmatch(deblank(M_.endo_names(k2,:)),deblank(options_.varobs),'exact');
+ estim_params_.corrn_observable_correspondence(i,:)=[obsi1,obsi2]; %save correspondence
+ end
end
if np
xparam1 = [xparam1; estim_params_.param_vals(:,2)];
--
GitLab