Merge branch 'MoM_optim_gradient' into 'master'

Method of Moments: use analytical Jacobian in gradient-based optimizers (GMM only)

See merge request Dynare/dynare!1801

matlab/method_of_moments/method_of_moments.m

@@ -93,6 +93,7 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_,
 % - [ ] 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
+% - [ ] add analytic_jacobian option for mode_compute 4 and 101
 % -------------------------------------------------------------------------
 % Step 0: Check if required structures and options exist
 % -------------------------------------------------------------------------
@@ -165,9 +166,10 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
         error('method_of_moments: perturbation orders higher than 3 are not implemented for GMM estimation, try using SMM.\n');
     end
     options_mom_.mom = set_default_option(options_mom_.mom,'analytic_standard_errors',false);       % compute standard errors numerically (0) or analytically (1). Analytical derivatives are only available for GMM.
+    options_mom_.mom = set_default_option(options_mom_.mom,'analytic_jacobian',false);              % use analytic Jacobian in optimization, only available for GMM and gradient-based optimizers    
 end
-options_mom_.mom.compute_derivs = false;% flag to compute derivs in objective function (needed for analytic standard errors with GMM)
+% initialize flag to compute derivs in objective function (needed for GMM with either analytic_standard_errors or analytic_jacobian )
-
+options_mom_.mom.compute_derivs = false;
     
 % General options that can be set by the user in the mod file, otherwise default values are provided
 options_mom_ = set_default_option(options_mom_,'dirname',M_.dname);    % specify directory in which to store estimation output [not yet working]
@@ -330,6 +332,10 @@ options_mom_.analytic_derivation_mode = 0; % needed by get_perturbation_params_d
 
 options_mom_.vector_output= false;           % specifies whether the objective function returns a vector
 
+optimizer_vec=[options_mom_.mode_compute;num2cell(options_mom_.additional_optimizer_steps)]; % at each stage one can possibly use different optimizers sequentially
+
+analytic_jacobian_optimizers = [1, 3, 13]; %these are currently supported, see to-do list
+
 % -------------------------------------------------------------------------
 % Step 1d: Other options that need to be initialized
 % -------------------------------------------------------------------------
@@ -708,6 +714,7 @@ end
 % Step 6: checks for objective function at initial parameters
 % -------------------------------------------------------------------------
 objective_function = str2func('method_of_moments_objective_function');
+
 try
     % Check for NaN or complex values of moment-distance-funtion evaluated
     % at initial parameters and identity weighting matrix    
@@ -757,7 +764,7 @@ end
 if options_mom_.mom.penalized_estimator
     fprintf('\n  - penalized estimation using deviation from prior mean and weighted with prior precision');
 end
-optimizer_vec=[options_mom_.mode_compute;num2cell(options_mom_.additional_optimizer_steps)]; % at each stage one can possibly use different optimizers sequentially
+
 for i = 1:length(optimizer_vec)
     if i == 1
         str = '- optimizer (mode_compute';
@@ -807,6 +814,9 @@ for i = 1:length(optimizer_vec)
     if options_mom_.silent_optimizer
         fprintf(' (silent)');
     end
+    if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_jacobian && ismember(optimizer_vec{i},analytic_jacobian_optimizers)
+        fprintf(' (using analytical Jacobian)');
+    end
 end
 fprintf('\n  - perturbation order:        %d', options_mom_.order)
 if options_mom_.order > 1 && options_mom_.pruning
@@ -828,6 +838,7 @@ if size(options_mom_.mom.weighting_matrix,1)>1 && ~(any(strcmpi('diagonal',optio
     fprintf('\nYou did not specify the use of an optimal or diagonal weighting matrix. There is no point in running an iterated method of moments.\n')
 end
 
+optim_opt0 = options_mom_.optim_opt; % store original options set by user
 for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
     fprintf('Estimation stage %u\n',stage_iter);
     Woptflag = false;
@@ -873,16 +884,36 @@ for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
     end
 
     for optim_iter= 1:length(optimizer_vec)
+        options_mom_.current_optimizer = optimizer_vec{optim_iter};
         if optimizer_vec{optim_iter}==0
             xparam1=xparam0; %no minimization, evaluate objective at current values
             fval = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_);
         else
             if optimizer_vec{optim_iter}==13
-                options_mom_.vector_output = true;
+                options_mom_.vector_output = true;                
             else
-                options_mom_.vector_output = false;
+                options_mom_.vector_output = false;                
             end
-            [xparam1, fval, exitflag] = dynare_minimize_objective(objective_function, xparam0, optimizer_vec{optim_iter}, options_mom_, [Bounds.lb Bounds.ub], bayestopt_laplace.name, bayestopt_laplace, [],...
+            if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_jacobian && ismember(optimizer_vec{optim_iter},analytic_jacobian_optimizers) %do this only for gradient-based optimizers
+                options_mom_.mom.compute_derivs = true;
+                objective_function1 = str2func('method_of_moments_objective_function_gradient_helper');
+                switch optimizer_vec{optim_iter}
+                    case 1
+                        options_mom_.optim_opt = [optim_opt0, ',''GradObj'',''on'''];   % make sure GradObj option is on for fmincon
+                    case 3
+                        options_mom_.optim_opt = [optim_opt0, ',''GradObj'',''on'''];   % make sure GradObj option is on for fmincon
+                    case 13
+                        options_mom_.optim_opt = [optim_opt0, ',''Jacobian'',''on'''];  % make sure Jacobian option is on for lsqnonlin
+                end
+                if strcmp(options_mom_.optim_opt(1),',')
+                    options_mom_.optim_opt(1) = []; %remove the comma if optim_opt was empty
+                end
+            else
+                options_mom_.mom.compute_derivs = false;
+                objective_function1 = objective_function;
+            end
+            
+            [xparam1, fval, exitflag] = dynare_minimize_objective(objective_function1, xparam0, optimizer_vec{optim_iter}, options_mom_, [Bounds.lb Bounds.ub], bayestopt_laplace.name, bayestopt_laplace, [],...
                                                                   Bounds, oo_, estim_params_, M_, options_mom_);
             if options_mom_.vector_output
                 fval = fval'*fval;

matlab/method_of_moments/method_of_moments_objective_function.m

@@ -1,4 +1,4 @@
-function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
+function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_, df] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
 % [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
 % -------------------------------------------------------------------------
 % This function evaluates the objective function for GMM/SMM estimation
@@ -21,9 +21,12 @@ function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_] = method_o
 %      - mom.model_moments       [numMom x 1] vector with model moments
 %      - mom.Q                   value of the quadratic form of the moment difference
 %   o M_:                       Matlab's structure describing the model
+%   o options_mom_:             structure information about all settings (specified by the user, preprocessor, and taken from global options_)
+%   o df:                       analytical parameter Jacobian of the quadratic form of the moment difference (for GMM only)
 % -------------------------------------------------------------------------
 % This function is called by
 %  o method_of_moments.m
+%  o dynare_minimize_objective.m
 % -------------------------------------------------------------------------
 % This function calls
 %  o check_bounds_and_definiteness_estimation
@@ -56,7 +59,15 @@ function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_] = method_o
 %------------------------------------------------------------------------------
 % 0. Initialization of the returned variables and others...
 %------------------------------------------------------------------------------
-
+if options_mom_.vector_output == 1
+    if options_mom_.mom.penalized_estimator
+        df           = nan(size(oo_.mom.data_moments,1)+length(xparam1),length(xparam1));
+    else
+        df           = nan(size(oo_.mom.data_moments,1),length(xparam1));
+    end
+else
+    df           = nan(1,length(xparam1));
+end
 junk1        = [];
 junk2        = [];
 
@@ -112,7 +123,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
     %--------------------------------------------------------------------------
     % 3. Set up pruned state-space system and compute model moments
     %--------------------------------------------------------------------------
-    if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors        
+    if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
         indpmodel = []; %initialize index for model parameters
         if ~isempty(estim_params_.param_vals)
             indpmodel = estim_params_.param_vals(:,1); %values correspond to parameters declaration order, row number corresponds to order in estimated_params
@@ -146,7 +157,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
         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
+        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;
@@ -156,12 +167,12 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
     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
+            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';
-            if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
+            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)';
@@ -170,7 +181,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
         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
+        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;
@@ -179,12 +190,12 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
     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
+            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
+            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)])...
@@ -194,7 +205,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
         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
+        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
@@ -265,6 +276,30 @@ else
     end
 end
 
+if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian
+    if options_mom_.mom.penalized_estimator
+        dxparam1 = eye(length(xparam1));
+    end
+    
+    for jp=1:length(xparam1)
+        dmoments_difference = - oo_.mom.model_moments_params_derivs(:,jp);
+        dresiduals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*oo_.mom.Sw*dmoments_difference;
+        
+        if options_mom_.vector_output == 1 % lsqnonlin requires vector output            
+            if options_mom_.mom.penalized_estimator                
+                df(:,jp)=[dresiduals;dxparam1(:,jp)./sqrt(diag(oo_.prior.variance))];
+            else
+                df(:,jp) = dresiduals;
+            end
+        else
+            df(:,jp) = dresiduals'*residuals + residuals'*dresiduals;
+            if options_mom_.mom.penalized_estimator
+                df(:,jp)=df(:,jp)+(dxparam1(:,jp))'/oo_.prior.variance*(xparam1-oo_.prior.mean)+(xparam1-oo_.prior.mean)'/oo_.prior.variance*(dxparam1(:,jp));
+            end
+        end
+    end
+end
+
 
 end%main function end
 

matlab/method_of_moments/method_of_moments_objective_function_gradient_helper.m

@@ -0,0 +1,55 @@
+function [out1, out2] = method_of_moments_objective_function_gradient_helper(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
+% [out1, out2] = method_of_moments_objective_function_gradient_helper(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
+% -------------------------------------------------------------------------
+% This helper function evaluates the objective function for GMM/SMM estimation and
+% outputs the function value fval at first and the gradient df at second place,
+% needed for gradient-based optimizers if analytic_jacobian option is set
+% =========================================================================
+% INPUTS
+%   o xparam1:                  current value of estimated parameters as returned by set_prior()
+%   o Bounds:                   structure containing parameter bounds
+%   o oo_:                      structure for results
+%   o estim_params_:            structure describing the estimated_parameters
+%   o M_                        structure describing the model
+%   o options_mom_:             structure information about all settings (specified by the user, preprocessor, and taken from global options_)
+% -------------------------------------------------------------------------
+% OUTPUTS: dependent on the optimizer calling this function, see below
+% -------------------------------------------------------------------------
+% This function calls
+%  o method_of_moments_objective_function.m
+% =========================================================================
+% Copyright (C) 2020-2021 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+% -------------------------------------------------------------------------
+% This function is called by
+%  o method_of_moments.m
+%  o dynare_minimize_objective.m
+% -------------------------------------------------------------------------
+% This function calls
+%  o method_of_moments_objective_function
+% -------------------------------------------------------------------------
+% Author(s): 
+% o Willi Mutschler (willi@mutschler.eu)
+% =========================================================================
+[fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_, df] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_);
+
+% mode_compute=1|3|13 require the gradient as second output argument
+out1=fval;
+out2=df;
+
+end%main function end
+

matlab/method_of_moments/method_of_moments_standard_errors.m

@@ -25,8 +25,7 @@ function [SE_values, Asympt_Var] = method_of_moments_standard_errors(xparam, obj
 % This function calls:
 %  o get_the_name
 %  o get_error_message
-%  o GMM_objective_function
+%  o method_of_moments_objective_function
-%  o SMM_objective_function.m
 %  o method_of_moments_optimal_weighting_matrix  
 % =========================================================================
 % Copyright (C) 2020-2021 Dynare Team

tests/Makefile.am

@@ -56,6 +56,7 @@ MODFILES = \
 	estimation/method_of_moments/RBC/RBC_MoM_SMM_ME.mod \
 	estimation/method_of_moments/RBC/RBC_MoM_prefilter.mod \
 	estimation/method_of_moments/RBC/RBC_MoM_optimizer.mod \
+	estimation/method_of_moments/RBC/RBC_MoM_GMM_gradient_optim.mod \
 	estimation/method_of_moments/AFVRR/AFVRR_M0.mod \
 	estimation/method_of_moments/AFVRR/AFVRR_MFB.mod \
 	estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod \

tests/estimation/method_of_moments/AFVRR/AFVRR_M0.mod

@@ -256,11 +256,13 @@ dev_Q            = AndreasenEtAl.Q@{orderApp} - oo_.mom.Q;
 dev_datamoments  = str2double(AndreasenEtAl.moments@{orderApp}(:,5)) - oo_.mom.data_moments;
 dev_modelmoments = str2double(AndreasenEtAl.moments@{orderApp}(:,6)) - oo_.mom.model_moments;
 
+if ~isoctave %there is no table command in Octave
 table([AndreasenEtAl.Q@{orderApp} ; str2double(AndreasenEtAl.moments@{orderApp}(:,5)) ; str2double(AndreasenEtAl.moments@{orderApp}(:,6))],...
       [oo_.mom.Q                  ; oo_.mom.data_moments                              ; oo_.mom.model_moments                            ],...
       [dev_Q                      ; dev_datamoments                                   ; dev_modelmoments                                 ],...
       'VariableNames', {'Andreasen et al', 'Dynare', 'dev'},...
       'RowNames', ['Q'; strcat('Data_', M_.matched_moments(:,4)); strcat('Model_', M_.matched_moments(:,4))])
+end
 
 if norm(dev_modelmoments)> 1e-4
     error('Something wrong in the computation of moments at order @{orderApp}')

tests/estimation/method_of_moments/AFVRR/AFVRR_MFB.mod

@@ -257,11 +257,13 @@ dev_Q            = AndreasenEtAl.Q@{orderApp} - oo_.mom.Q;
 dev_datamoments  = str2double(AndreasenEtAl.moments@{orderApp}(:,5)) - oo_.mom.data_moments;
 dev_modelmoments = str2double(AndreasenEtAl.moments@{orderApp}(:,6)) - oo_.mom.model_moments;
 
+if ~isoctave %there is no table command in Octave
 table([AndreasenEtAl.Q@{orderApp} ; str2double(AndreasenEtAl.moments@{orderApp}(:,5)) ; str2double(AndreasenEtAl.moments@{orderApp}(:,6))],...
       [oo_.mom.Q                  ; oo_.mom.data_moments                              ; oo_.mom.model_moments                            ],...
       [dev_Q                      ; dev_datamoments                                   ; dev_modelmoments                                 ],...
       'VariableNames', {'Andreasen et al', 'Dynare', 'dev'},...
       'RowNames', ['Q'; strcat('Data_', M_.matched_moments(:,4)); strcat('Model_', M_.matched_moments(:,4))])
+end
 
 if norm(dev_modelmoments)> 1e-4
     warning('Something wrong in the computation of moments at order @{orderApp}')

tests/estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod

@@ -256,11 +256,13 @@ dev_Q            = AndreasenEtAl.Q@{orderApp} - oo_.mom.Q;
 dev_datamoments  = str2double(AndreasenEtAl.moments@{orderApp}(:,5)) - oo_.mom.data_moments;
 dev_modelmoments = str2double(AndreasenEtAl.moments@{orderApp}(:,6)) - oo_.mom.model_moments;
 
+if ~isoctave %there is no table command in Octave
 table([AndreasenEtAl.Q@{orderApp} ; str2double(AndreasenEtAl.moments@{orderApp}(:,5)) ; str2double(AndreasenEtAl.moments@{orderApp}(:,6))],...
       [oo_.mom.Q                  ; oo_.mom.data_moments                              ; oo_.mom.model_moments                            ],...
       [dev_Q                      ; dev_datamoments                                   ; dev_modelmoments                                 ],...
       'VariableNames', {'Andreasen et al', 'Dynare', 'dev'},...
       'RowNames', ['Q'; strcat('Data_', M_.matched_moments(:,4)); strcat('Model_', M_.matched_moments(:,4))])
+end
 
 if norm(dev_modelmoments)> 1e-4
     warning('Something wrong in the computation of moments at order @{orderApp}')

tests/estimation/method_of_moments/RBC/RBC_MoM_GMM_gradient_optim.mod

@@ -0,0 +1,125 @@
+% Test whether gradient-based optimizers are able to use analytical
+% Jacobian of moments in GMM estimation
+%
+% Copyright (C) 2021 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+% =========================================================================
+@#include "RBC_MoM_common.inc"
+
+shocks;
+var u_a; stderr 0.0072;        
+end;
+
+varobs c iv n;
+
+%--------------------------------------------------------------------------
+% Method of Moments Estimation
+%--------------------------------------------------------------------------
+matched_moments;
+c;
+n;
+iv;
+c*c;
+c*iv;
+iv*n;
+iv*iv;
+n*c;
+n*n;
+c*c(-1);
+n*n(-1);
+iv*iv(-1);
+end;
+
+
+@#for estimParams in [0, 1, 2]
+  clear estim_params_;
+
+  @#if estimParams == 0
+    estimated_params;
+        %DELTA,         0.025;
+        %BETTA,         0.984;
+        %B,             0.5;
+        %ETAc,          2;
+        ALFA,          0.667;
+        RHOA,          0.979;
+        stderr u_a,    0.0072;
+    end;
+  @#endif
+
+  @#if estimParams == 1
+    estimated_params;
+        %DELTA,         ,        0,           1;
+        %BETTA,         ,        0,           1;
+        %B,             ,        0,           1;
+        %ETAc,          ,        0,           10;
+        ALFA,          ,        0,           1;
+        RHOA,          ,        0,           1;
+        stderr u_a,    ,        0,           1;
+    end;
+  @#endif
+
+  @#if estimParams == 2
+    estimated_params;
+        %DELTA,         0.025,         0,           1,  normal_pdf, 0.02, 0.5;
+        %BETTA,         0.98,         0,           1,  beta_pdf, 0.90, 0.25;
+        %B,             0.45,         0,           1,  normal_pdf, 0.40, 0.5;
+        %ETAl,          1,            0,           10, normal_pdf, 0.25, 0.0.1;
+        %ETAc,          1.8,         0,           10, normal_pdf, 1.80, 0.5;
+        ALFA,          0.65,         0,           1,  normal_pdf, 0.60, 0.5;
+        RHOA,          0.95,         0,           1,  normal_pdf, 0.90, 0.5;
+        stderr u_a,    0.01,         0,           1,  normal_pdf, 0.01, 0.5;
+        %THETA,         3.48,          0,           10, normal_pdf, 0.25, 0.0.1;
+    end;
+  @#endif
+
+  estimated_params_init(use_calibration);
+  end;
+  
+  @#for optimizer in [1, 3, 13]
+    @#if estimParams == 2 && optimizer == 13
+        %skip due to buggy behavior in Octave
+        if ~isoctave
+    @#endif
+  
+  method_of_moments(
+          mom_method = GMM         % method of moments method; possible values: GMM|SMM
+        , datafile   = 'RBC_Andreasen_Data_2.mat' % name of filename with data
+        , order = 2                 % order of Taylor approximation in perturbation
+        , weighting_matrix = ['OPTIMAL']      % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename. Size of cell determines stages in iterated estimation, e.g. two state with ['DIAGONAL','OPTIMAL']
+        , nodisplay
+        , nograph
+        , mode_compute = @{optimizer}        % specifies the optimizer for minimization of moments distance
+        %, additional_optimizer_steps = [1 3 13]
+%        , optim = ('TolFun', 1e-6
+%                    ,'TolX', 1e-6
+%                    ,'MaxIter', 3000
+%                    ,'MaxFunEvals', 1D6
+%                    ,'UseParallel' , 1
+%                    ,'Jacobian' , 'on'                   
+%                    ,'GradObj','on'
+%                   )    % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
+        , silent_optimizer                  % run minimization of moments distance silently without displaying results or saving files in between        
+        , analytic_jacobian
+  );
+  
+    @#if estimParams == 2 && optimizer == 13
+        %skip due to buggy behavior in Octave
+        end
+    @#endif
+  @#endfor  
+  
+@#endfor

tests/estimation/method_of_moments/RBC/RBC_MoM_optimizer.mod

@@ -118,6 +118,12 @@ options_.solveopt.TolXConstraint=1e-3;
     end;
 
     @#for optimizer in OPTIMIZERS
+        
+        @#if estimParams == 2 && optimizer == 13
+            %skip due to buggy behavior in Octave
+        if ~isoctave
+        @#endif
+        
     method_of_moments(
           mom_method = GMM         % method of moments method; possible values: GMM|SMM
         , datafile   = 'RBC_Andreasen_Data_2.mat' % name of filename with data
@@ -139,6 +145,11 @@ options_.solveopt.TolXConstraint=1e-3;
         @#endif
         %, silent_optimizer                  % run minimization of moments distance silently without displaying results or saving files in between
     );
+    
+        @#if estimParams == 2 && optimizer == 13
+            %skip due to buggy behavior in Octave
+        end
+        @#endif
     @#endfor
 @#endfor