diff --git a/matlab/method_of_moments/method_of_moments.m b/matlab/method_of_moments/method_of_moments.m
index e4f3a3ffb2992efb3c077428ec640088d3d02bc8..eb5579f43a694597103effa5999d678cee2270df 100644
--- a/matlab/method_of_moments/method_of_moments.m
+++ b/matlab/method_of_moments/method_of_moments.m
@@ -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
diff --git a/matlab/method_of_moments/method_of_moments_objective_function.m b/matlab/method_of_moments/method_of_moments_objective_function.m
index cb069ef1820b40aba1aad04a8c778f8ec9fffa5f..99b7d91aec04d9ae591d4f907a16e9489f8b532c 100644
--- a/matlab/method_of_moments/method_of_moments_objective_function.m
+++ b/matlab/method_of_moments/method_of_moments_objective_function.m
@@ -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