diff --git a/matlab/method_of_moments/method_of_moments.m b/matlab/method_of_moments/method_of_moments.m
index b931bb0573b48f5f961195649d5f4d59272032e6..5ea218c2939230e4a6bd7aeaeab9879defe0697c 100644
--- a/matlab/method_of_moments/method_of_moments.m
+++ b/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;
+ end
+ 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_.vector_output = false;
+ options_mom_.mom.compute_derivs = false;
+ objective_function1 = objective_function;
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, [],...
+
+ [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;
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 41c3dcc79ccbd9192753aacf359d1356213b8234..3a8548e76105de126ddeb8fdb2980fc20cccd2bf 100644
--- a/matlab/method_of_moments/method_of_moments_objective_function.m
+++ b/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
diff --git a/matlab/method_of_moments/method_of_moments_objective_function_gradient_helper.m b/matlab/method_of_moments/method_of_moments_objective_function_gradient_helper.m
new file mode 100644
index 0000000000000000000000000000000000000000..f1c88af499ae2d8367c2a62ddd9a65b3a96230cb
--- /dev/null
+++ b/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
+
diff --git a/matlab/method_of_moments/method_of_moments_standard_errors.m b/matlab/method_of_moments/method_of_moments_standard_errors.m
index 459f583f7aade842d0804353471da31f594c14ea..ea302d52d57d7252e663ee4c18787844fc8daa80 100644
--- a/matlab/method_of_moments/method_of_moments_standard_errors.m
+++ b/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 SMM_objective_function.m
+% o method_of_moments_objective_function
% o method_of_moments_optimal_weighting_matrix
% =========================================================================
% Copyright (C) 2020-2021 Dynare Team
diff --git a/tests/Makefile.am b/tests/Makefile.am
index ef844e1bdc90b554daf1793c1585feb31669af1b..37825c1cced223b2399b7315bdd6f6876fb4adfe 100644
--- a/tests/Makefile.am
+++ b/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 \
diff --git a/tests/estimation/method_of_moments/AFVRR/AFVRR_M0.mod b/tests/estimation/method_of_moments/AFVRR/AFVRR_M0.mod
index 8e51ac51368237c54d02210941407fdfd3a8db08..cd19168cc962c612acd082a692646cf118fec7f8 100644
--- a/tests/estimation/method_of_moments/AFVRR/AFVRR_M0.mod
+++ b/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}')
diff --git a/tests/estimation/method_of_moments/AFVRR/AFVRR_MFB.mod b/tests/estimation/method_of_moments/AFVRR/AFVRR_MFB.mod
index 450739ad3bd153855806d8bcbab513d554eb33b7..535dc847016bd2b7c753e8b8c3e067b5b9b8eefc 100644
--- a/tests/estimation/method_of_moments/AFVRR/AFVRR_MFB.mod
+++ b/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}')
diff --git a/tests/estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod b/tests/estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod
index 9c069d3a3dff9dfb87ce7e40f2678f27ce1364aa..d86c7b35374de54866af804af4f4bfaa3a10e32f 100644
--- a/tests/estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod
+++ b/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}')
diff --git a/tests/estimation/method_of_moments/RBC/RBC_MoM_GMM_gradient_optim.mod b/tests/estimation/method_of_moments/RBC/RBC_MoM_GMM_gradient_optim.mod
new file mode 100644
index 0000000000000000000000000000000000000000..39a58b6740efad0b5c656ddff9051d3a1ceea316
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/RBC/RBC_MoM_optimizer.mod b/tests/estimation/method_of_moments/RBC/RBC_MoM_optimizer.mod
index 08b1ed882fabf4ec1bea83047d394e736400ca63..b9d9dc96c612671f4735c677e62be67453fcacd3 100644
--- a/tests/estimation/method_of_moments/RBC/RBC_MoM_optimizer.mod
+++ b/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