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 <willi@mutschler.eu>
(cherry picked from commit fdba1170)

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
-
-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
+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 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) ,'})$'];
-        end
-    end
-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);
+% check for duplicate moment conditions
+for jm=1:size(M_.matched_moments,1)
+    % 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
+
+% 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
+
+% 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: %s.\n',num2str(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
-%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');
+    % 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);
+    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

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;
+    for jm = 1:size(M_.matched_moments,1)
         % 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_.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 )
-            oo_.mom.model_moments_params_derivs(offset+1:E_y_nbr,:) = pruned_state_space.dE_y(options_mom_.mom.index.E_y,:);
+                oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dE_y(idx1,:);
             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 (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
-            E_yy = pruned_state_space.Var_y;
+                    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 )
-                dE_yy = pruned_state_space.dVar_y;
+                        oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dVar_y(idx1,idx2,:);
                     end
                 else
-            E_yy = pruned_state_space.Var_y + pruned_state_space.E_y*pruned_state_space.E_y';
+                    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 )
-                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)';
+                            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            
-        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
+            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
-            E_yyt = pruned_state_space.Var_yi;
+                    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 )
-                dE_yyt = pruned_state_space.dVar_yi;
+                        oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dVar_yi(idx1,idx2,lag,:);
                     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)]);
+                    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 )
-                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)]);
+                            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
-        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
     
+    
 elseif strcmp(options_mom_.mom.mom_method,'SMM')
     %------------------------------------------------------------------------------
     % 3. Compute Moments of the model solution for normal innovations