Verified Commit fdba1170 authored by Willi Mutschler's avatar Willi Mutschler Committed by Stéphane Adjemian
Browse files

MoM: Refactor check for duplicate moments, fix display of results



As we have a working interface now, this commit improves the provisional handling of finding duplicate moments. Previously, indices for GMM were created, but this is not really needed. This commit cleans this up and similar to SMM makes use of the matched_moments block.

As a by-product of the previous provisonal handling higher-order moments for SMM where not correctly displayed as no labels were created. This is now fixed. The comparison of data moments and estimated model moments is also in the same ordering as the inputed orthogonality conditions in matched_moments.
Signed-off-by: Willi Mutschler's avatarWilli Mutschler <willi@mutschler.eu>
parent 4d590f8a
......@@ -91,7 +91,6 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_,
% - [ ] add option to use autocorrelations (we have useautocorr in identification toolbox already)
% - [ ] SMM with extended path
% - [ ] deal with measurement errors (once @wmutschl has implemented this in identification toolbox)
% - [ ] improve check for duplicate moments by using the cellfun and unique functions
% - [ ] dirname option to save output to different directory not yet implemented
% The TeX option crashes MATLAB R2014a run with "-nodisplay" option
......@@ -384,85 +383,66 @@ for i=1:options_mom_.obs_nbr
end
% -------------------------------------------------------------------------
% Step 2: Checks and transformations for matched moments structure (preliminary)
% Step 2: Checks and transformations for matched_moments structure
% -------------------------------------------------------------------------
M_.matched_moments_orig = M_.matched_moments; %save original structure
% Initialize indices
options_mom_.mom.index.E_y = false(options_mom_.obs_nbr,1); %unconditional first order product moments
options_mom_.mom.index.E_yy = false(options_mom_.obs_nbr,options_mom_.obs_nbr); %unconditional second order product moments
options_mom_.mom.index.E_yyt = false(options_mom_.obs_nbr,options_mom_.obs_nbr,0); %unconditional temporal second order product moments
options_mom_.mom.index.E_y_pos = zeros(options_mom_.obs_nbr,1); %position in matched moments block
options_mom_.mom.index.E_yy_pos = zeros(options_mom_.obs_nbr,options_mom_.obs_nbr); %position in matched moments block
options_mom_.mom.index.E_yyt_pos = zeros(options_mom_.obs_nbr,options_mom_.obs_nbr,0); %position in matched moments block
% higher-order product moments not supported yet for GMM
if strcmp(options_mom_.mom.mom_method, 'GMM') && any(cellfun(@sum,M_.matched_moments(:,3))> 2)
error('method_of_moments: GMM does not yet support product moments higher than 2. Change row %d in ''matched_moments'' block.',jm);
end
% check for duplicate moment conditions
for jm=1:size(M_.matched_moments,1)
% higher-order product moments not supported yet for GMM
if strcmp(options_mom_.mom.mom_method, 'GMM') && sum(M_.matched_moments{jm,3}) > 2
error('method_of_moments: GMM does not yet support product moments higher than 2. Change row %d in ''matched_moments'' block.',jm);
end
% Check if declared variables are also observed (needed as otherwise the dataset variables won't coincide)
if any(~ismember(oo_.dr.inv_order_var(M_.matched_moments{jm,1})', oo_.dr.obs_var))
error('method_of_moments: Variables in row %d in ''matched_moments'' block need to be declared as VAROBS.', jm)
end
if strcmp(options_mom_.mom.mom_method, 'GMM')
% Check (for now) that only lags are declared
if any(M_.matched_moments{jm,2}>0)
error('method_of_moments: Leads in row %d in the ''matched_moments'' block are not supported for GMM, shift the moments and declare only lags.', jm)
end
% Check (for now) that first declared variable has zero lag
if M_.matched_moments{jm,2}(1)~=0
error('method_of_moments: The first variable declared in row %d in the ''matched_moments'' block is not allowed to have a lead or lag for GMM;\n reorder the variables in the row such that the first variable has zero lag!',jm)
end
end
vars = oo_.dr.inv_order_var(M_.matched_moments{jm,1})';
if sum(M_.matched_moments{jm,3}) == 1
% First-order product moment
vpos = (oo_.dr.obs_var == vars);
options_mom_.mom.index.E_y(vpos,1) = true;
options_mom_.mom.index.E_y_pos(vpos,1) = jm;
M_.matched_moments{jm,4}=['E(',M_.endo_names{M_.matched_moments{jm,1}},')'];
M_.matched_moments{jm,5}=['$E(',M_.endo_names_tex{M_.matched_moments{jm,1}},')$'];
elseif sum(M_.matched_moments{jm,3}) == 2
% Second-order product moment
idx1 = (oo_.dr.obs_var == vars(1));
idx2 = (oo_.dr.obs_var == vars(2));
lag1 = M_.matched_moments{jm,2}(1);
lag2 = M_.matched_moments{jm,2}(2);
if lag1==0 && lag2==0 % contemporaneous covariance matrix
options_mom_.mom.index.E_yy(idx1,idx2) = true;
options_mom_.mom.index.E_yy(idx2,idx1) = true;
options_mom_.mom.index.E_yy_pos(idx1,idx2) = jm;
options_mom_.mom.index.E_yy_pos(idx2,idx1) = jm;
M_.matched_moments{jm,4}=['E(',M_.endo_names{M_.matched_moments{jm,1}(1)},',',M_.endo_names{M_.matched_moments{jm,1}(2)},')'];
M_.matched_moments{jm,5}=['$E({',M_.endo_names_tex{M_.matched_moments{jm,1}(1)},'}_t,{',M_.endo_names_tex{M_.matched_moments{jm,1}(1)},'}_t)$'];
elseif lag1==0 && lag2 < 0
options_mom_.mom.index.E_yyt(idx1,idx2,-lag2) = true;
options_mom_.mom.index.E_yyt_pos(idx1,idx2,-lag2) = jm;
M_.matched_moments{jm,4}=['E(',M_.endo_names{M_.matched_moments{jm,1}(1)},',',M_.endo_names{M_.matched_moments{jm,1}(2)},'(',num2str(lag2),'))'];
M_.matched_moments{jm,5}=['$E({',M_.endo_names_tex{M_.matched_moments{jm,1}(1)},'}_t\times{',M_.endo_names_tex{M_.matched_moments{jm,1}(1)},'_{t',num2str(lag2) ,'})$'];
% expand powers to vector of ones
if any(M_.matched_moments{jm,3}>1)
tmp1=[]; tmp2=[]; tmp3=[];
for jjm=1:length(M_.matched_moments{jm,3})
tmp1 = [tmp1 repmat(M_.matched_moments{jm,1}(jjm),[1 M_.matched_moments{jm,3}(jjm)]) ];
tmp2 = [tmp2 repmat(M_.matched_moments{jm,2}(jjm),[1 M_.matched_moments{jm,3}(jjm)]) ];
tmp3 = [tmp3 repmat(1,[1 M_.matched_moments{jm,3}(jjm)]) ];
end
M_.matched_moments{jm,1} = tmp1;
M_.matched_moments{jm,2} = tmp2;
M_.matched_moments{jm,3} = tmp3;
end
% shift time structure to focus only on lags
M_.matched_moments{jm,2} = M_.matched_moments{jm,2} - max(M_.matched_moments{jm,2});
% sort such that t=0 variable comes first
[M_.matched_moments{jm,2},idx_sort] = sort(M_.matched_moments{jm,2},'descend');
M_.matched_moments{jm,1} = M_.matched_moments{jm,1}(idx_sort);
M_.matched_moments{jm,3} = M_.matched_moments{jm,3}(idx_sort);
end
%Remove duplicate elements
UniqueMomIdx = [nonzeros(options_mom_.mom.index.E_y_pos); nonzeros(tril(options_mom_.mom.index.E_yy_pos)); nonzeros(options_mom_.mom.index.E_yyt_pos)];
DuplicateMoms = setdiff(1:size(M_.matched_moments,1),UniqueMomIdx);
if ~isempty(DuplicateMoms)
fprintf('Found and removed duplicate declared moments in ''matched_moments'' block in rows: %s.\n',num2str(DuplicateMoms))
% find duplicate rows in cell array by making groups according to powers as we can then use cell2mat for the unique function
powers = cellfun(@sum,M_.matched_moments(:,3))';
UniqueMomIdx = [];
for jpow = unique(powers)
idx1 = find(powers==jpow);
[~,idx2] = unique(cell2mat(M_.matched_moments(idx1,:)),'rows');
UniqueMomIdx = [UniqueMomIdx idx1(idx2)];
end
%reorder M_.matched_moments to be compatible with options_mom_.mom.index
M_.matched_moments = M_.matched_moments(UniqueMomIdx,:);
if strcmp(options_mom_.mom.mom_method,'SMM')
options_mom_.mom=rmfield(options_mom_.mom,'index');
% remove duplicate elements
DuplicateMoms = setdiff(1:size(M_.matched_moments_orig,1),UniqueMomIdx);
if ~isempty(DuplicateMoms)
fprintf('Found and removed duplicate declared moments in ''matched_moments'' block in rows:\n %s.\n',num2str(DuplicateMoms))
fprintf('Dynare will continue with remaining moment conditions\n');
end
if strcmp(options_mom_.mom.mom_method, 'SMM')
% for SMM we can keep the original structure but get rid of duplicate moments
M_.matched_moments = M_.matched_moments_orig(sort(UniqueMomIdx),:);
elseif strcmp(options_mom_.mom.mom_method, 'GMM')
% for GMM we use the transformed matched_moments structure
M_.matched_moments = M_.matched_moments(sort(UniqueMomIdx),:);
end
% Check if both prefilter and first moments were specified
options_mom_.mom.first_moment_indicator = find(cellfun(@(x) sum(abs(x))==1,M_.matched_moments(:,3)))';
if options_mom_.prefilter && ~isempty(options_mom_.mom.first_moment_indicator)
fprintf('Centered moments requested (prefilter option is set); therefore, ignore declared first moments in ''matched_moments'' block in rows: %u.\n',options_mom_.mom.first_moment_indicator');
M_.matched_moments(options_mom_.mom.first_moment_indicator,:)=[]; %remove first moments entries
options_mom_.mom.first_moment_indicator = [];
first_moment_indicator = find(cellfun(@(x) sum(abs(x))==1,M_.matched_moments(:,3)))';
if options_mom_.prefilter && ~isempty(first_moment_indicator)
fprintf('Centered moments requested (prefilter option is set); therefore, ignore declared first moments in ''matched_moments'' block.\n');
M_.matched_moments(first_moment_indicator,:)=[]; %remove first moments entries
end
options_mom_.mom.mom_nbr = size(M_.matched_moments,1);
......@@ -972,15 +952,39 @@ end
% -------------------------------------------------------------------------
% Step 9: Display estimation results
% -------------------------------------------------------------------------
title = ['Data moments and model moments (',options_mom_.mom.mom_method,')'];
headers = {'Moment','Data','Model','% dev. target'};
labels= M_.matched_moments(:,4);
data_mat=[oo_.mom.data_moments oo_.mom.model_moments 100*abs((oo_.mom.model_moments-oo_.mom.data_moments)./oo_.mom.data_moments)];
title = ['Comparison of data moments and model moments (',options_mom_.mom.mom_method,')'];
headers = {'Moment','Data','Model'};
for jm = 1:size(M_.matched_moments,1)
lables_tmp = 'E[';
lables_tmp_tex = 'E \left[ ';
for jvar = 1:length(M_.matched_moments{jm,1})
lables_tmp = [lables_tmp M_.endo_names{M_.matched_moments{jm,1}(jvar)}];
lables_tmp_tex = [lables_tmp_tex, '{', M_.endo_names_tex{M_.matched_moments{jm,1}(jvar)}, '}'];
if M_.matched_moments{jm,2}(jvar) ~= 0
lables_tmp = [lables_tmp, '(', num2str(M_.matched_moments{jm,2}(jvar)), ')'];
lables_tmp_tex = [lables_tmp_tex, '_{t', num2str(M_.matched_moments{jm,2}(jvar)), '}'];
else
lables_tmp_tex = [lables_tmp_tex, '_{t}'];
end
if M_.matched_moments{jm,3}(jvar) > 1
lables_tmp = [lables_tmp, '^', num2str(M_.matched_moments{jm,3}(jvar))];
lables_tmp_tex = [lables_tmp_tex, '^{', num2str(M_.matched_moments{jm,3}(jvar)) '}'];
end
if jvar == length(M_.matched_moments{jm,1})
lables_tmp = [lables_tmp, ']'];
lables_tmp_tex = [lables_tmp_tex, ' \right]'];
else
lables_tmp = [lables_tmp, '*'];
lables_tmp_tex = [lables_tmp_tex, ' \times '];
end
end
labels{jm,1} = lables_tmp;
labels_TeX{jm,1} = lables_tmp_tex;
end
data_mat=[oo_.mom.data_moments oo_.mom.model_moments ];
dyntable(options_mom_, title, headers, labels, data_mat, cellofchararraymaxlength(labels)+2, 10, 7);
if options_mom_.TeX
lh = cellofchararraymaxlength(labels)+2;
labels_TeX = M_.matched_moments(:,5);
dyn_latex_table(M_, options_mom_, title, 'sim_corr_matrix', headers, labels_TeX, data_mat, lh, 10, 7);
if options_mom_.TeX
dyn_latex_table(M_, options_mom_, title, ['comparison_moments_', options_mom_.mom.mom_method], headers, labels_TeX, data_mat, cellofchararraymaxlength(labels)+2, 10, 7);
end
if options_mom_.mode_check.status
......
......@@ -154,65 +154,55 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
end
oo_.mom.model_moments = NaN(options_mom_.mom.mom_nbr,1);
offset = 0;
% First moments
if ~options_mom_.prefilter && isfield(options_mom_.mom.index,'E_y') && nnz(options_mom_.mom.index.E_y) > 0
E_y = pruned_state_space.E_y;
E_y_nbr = nnz(options_mom_.mom.index.E_y);
oo_.mom.model_moments(offset+1:E_y_nbr,1) = E_y(options_mom_.mom.index.E_y);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
oo_.mom.model_moments_params_derivs(offset+1:E_y_nbr,:) = pruned_state_space.dE_y(options_mom_.mom.index.E_y,:);
end
offset = offset + E_y_nbr;
end
% Second moments
% Contemporaneous covariance
if isfield(options_mom_.mom.index,'E_yy') && nnz(options_mom_.mom.index.E_yy) > 0
if options_mom_.prefilter
E_yy = pruned_state_space.Var_y;
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
dE_yy = pruned_state_space.dVar_y;
end
else
E_yy = pruned_state_space.Var_y + pruned_state_space.E_y*pruned_state_space.E_y';
for jm = 1:size(M_.matched_moments,1)
% First moments
if ~options_mom_.prefilter && (sum(M_.matched_moments{jm,3}) == 1)
idx1 = (oo_.dr.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}) );
oo_.mom.model_moments(jm,1) = pruned_state_space.E_y(idx1);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
dE_yy = pruned_state_space.dVar_y;
for jp=1:totparam_nbr
dE_yy(:,:,jp) = dE_yy(:,:,jp) + pruned_state_space.dE_y(:,jp)*pruned_state_space.E_y' + pruned_state_space.E_y*pruned_state_space.dE_y(:,jp)';
end
oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dE_y(idx1,:);
end
end
E_yy_nbr = nnz(tril(options_mom_.mom.index.E_yy));
oo_.mom.model_moments(offset+(1:E_yy_nbr),1) = E_yy(tril(options_mom_.mom.index.E_yy));
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
oo_.mom.model_moments_params_derivs(offset+(1:E_yy_nbr),:) = reshape(dE_yy(repmat(tril(options_mom_.mom.index.E_yy),[1 1 totparam_nbr])),E_yy_nbr,totparam_nbr);
end
offset = offset + E_yy_nbr;
end
% Lead/lags covariance
if isfield(options_mom_.mom.index,'E_yyt') && nnz(options_mom_.mom.index.E_yyt) > 0
if options_mom_.prefilter
E_yyt = pruned_state_space.Var_yi;
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
dE_yyt = pruned_state_space.dVar_yi;
end
else
E_yyt = pruned_state_space.Var_yi + repmat(pruned_state_space.E_y*pruned_state_space.E_y',[1 1 size(pruned_state_space.Var_yi,3)]);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
dE_yyt = pruned_state_space.dVar_yi;
for jp=1:totparam_nbr
dE_yyt(:,:,:,jp) = dE_yyt(:,:,:,jp) + repmat(pruned_state_space.dE_y(:,jp)*pruned_state_space.E_y',[1 1 size(pruned_state_space.Var_yi,3)])...
+ repmat(pruned_state_space.E_y*pruned_state_space.dE_y(:,jp)',[1 1 size(pruned_state_space.Var_yi,3)]);
end
% Second moments
if (sum(M_.matched_moments{jm,3}) == 2)
idx1 = (oo_.dr.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(1)) );
idx2 = (oo_.dr.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(2)) );
if nnz(M_.matched_moments{jm,2}) == 0
% Covariance
if options_mom_.prefilter
oo_.mom.model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dVar_y(idx1,idx2,:);
end
else
oo_.mom.model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2) + pruned_state_space.E_y(idx1)*pruned_state_space.E_y(idx2)';
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
for jp=1:totparam_nbr
oo_.mom.model_moments_params_derivs(jm,jp) = pruned_state_space.dVar_y(idx1,idx2,jp) + pruned_state_space.dE_y(idx1,jp)*pruned_state_space.E_y(idx2)' + pruned_state_space.E_y(idx1)*pruned_state_space.dE_y(idx2,jp)';
end
end
end
else
% Autocovariance
lag = -M_.matched_moments{jm,2}(2); %note that leads/lags in matched_moments are transformed such that first entry is always 0 and the second is a lag
if options_mom_.prefilter
oo_.mom.model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dVar_yi(idx1,idx2,lag,:);
end
else
oo_.mom.model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag) + pruned_state_space.E_y(idx1)*pruned_state_space.E_y(idx2)';
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
for jp=1:totparam_nbr
oo_.mom.model_moments_params_derivs(jm,jp) = vec( pruned_state_space.dVar_yi(idx1,idx2,lag,jp) + pruned_state_space.dE_y(idx1,jp)*pruned_state_space.E_y(idx2)' + pruned_state_space.E_y(idx1)*pruned_state_space.dE_y(idx2,jp)');
end
end
end
end
end
E_yyt_nbr = nnz(options_mom_.mom.index.E_yyt);
oo_.mom.model_moments(offset+(1:E_yyt_nbr),1) = E_yyt(options_mom_.mom.index.E_yyt);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
oo_.mom.model_moments_params_derivs(offset+(1:E_yyt_nbr),:) = reshape(dE_yyt(repmat(options_mom_.mom.index.E_yyt,[1 1 1 totparam_nbr])),E_yyt_nbr,totparam_nbr);
end
end
end
elseif strcmp(options_mom_.mom.mom_method,'SMM')
%------------------------------------------------------------------------------
% 3. Compute Moments of the model solution for normal innovations
......
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