Commit 7518072e authored by Johannes Pfeifer's avatar Johannes Pfeifer
Browse files

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
-
parent ebc7d6f6
......@@ -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
......
......@@ -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;'])
......
......@@ -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;
......@@ -554,7 +572,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), ...
......@@ -653,7 +671,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);']);
......@@ -799,7 +817,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,:)), ...
......@@ -1100,7 +1118,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)
......@@ -1131,7 +1150,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');
......@@ -1168,7 +1187,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
......
......@@ -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]);
......
......@@ -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);
......
......@@ -422,6 +422,7 @@ oo_.exo_det_simul = [];
M_.params = [];
M_.endo_histval = [];
M_.Correlation_matrix = [];
M_.Correlation_matrix_ME = [];
% homotopy
options_.homotopy_mode = 0;
......
......@@ -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
......
......@@ -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,:)))...
......
......@@ -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
......@@ -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)];
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment