diff --git a/matlab/method_of_moments/method_of_moments.m b/matlab/method_of_moments/method_of_moments.m
index 5d86fb93d..9d6736116 100644
--- a/matlab/method_of_moments/method_of_moments.m
+++ b/matlab/method_of_moments/method_of_moments.m
@@ -63,7 +63,7 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_,
% o set_all_parameters.m
% o test_for_deep_parameters_calibration.m
% =========================================================================
-% Copyright (C) 2020 Dynare Team
+% Copyright (C) 2020-2021 Dynare Team
%
% This file is part of Dynare.
%
@@ -92,6 +92,7 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_,
% - [ ] 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
% -------------------------------------------------------------------------
% Step 0: Check if required structures and options exist
% -------------------------------------------------------------------------
@@ -133,9 +134,9 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
options_mom_.mom = set_default_option(options_mom_.mom,'bartlett_kernel_lag',20); % bandwith in optimal weighting matrix
options_mom_.mom = set_default_option(options_mom_.mom,'penalized_estimator',false); % include deviation from prior mean as additional moment restriction and use prior precision as weight
options_mom_.mom = set_default_option(options_mom_.mom,'verbose',false); % display and store intermediate estimation results
- options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix',{'DIAGONAL'; 'DIAGONAL'}); % weighting matrix in moments distance objective function at each iteration of estimation; cell of strings with
+ options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix',{'DIAGONAL'; 'DIAGONAL'}); % weighting matrix in moments distance objective function at each iteration of estimation;
% possible values are 'OPTIMAL', 'IDENTITY_MATRIX' ,'DIAGONAL' or a filename. Size of cell determines stages in iterated estimation.
- options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix_scaling_factor',1); % scaling of weighting matrix
+ options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix_scaling_factor',1); % scaling of weighting matrix in objective function
options_mom_.mom = set_default_option(options_mom_.mom,'se_tolx',1e-5); % step size for numerical computation of standard errors
options_mom_ = set_default_option(options_mom_,'order',1); % order of Taylor approximation in perturbation
options_mom_ = set_default_option(options_mom_,'pruning',false); % use pruned state space system at higher-order
@@ -169,6 +170,7 @@ options_mom_.mom.compute_derivs = false;% flag to compute derivs in objective fu
% 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]
options_mom_ = set_default_option(options_mom_,'graph_format','eps'); % specify the file format(s) for graphs saved to disk
options_mom_ = set_default_option(options_mom_,'nodisplay',false); % do not display the graphs, but still save them to disk
options_mom_ = set_default_option(options_mom_,'nograph',false); % do not create graphs (which implies that they are not saved to the disk nor displayed)
@@ -200,7 +202,7 @@ options_mom_ = set_default_option(options_mom_,'optim_opt',[]);
options_mom_ = set_default_option(options_mom_,'silent_optimizer',false); % run minimization of moments distance silently without displaying results or saving files in between
% Check plot options that can be set by the user in the mod file, otherwise default values are provided
options_mom_.mode_check.nolik = false; % we don't do likelihood (also this initializes mode_check substructure)
-options_mom_.mode_check = set_default_option(options_mom_.mode_check,'status',false); % plot the target function for values around the computed minimum for each estimated parameter in turn. This is helpful to diagnose problems with the optimizer.
+options_mom_.mode_check = set_default_option(options_mom_.mode_check,'status',false); % plot the target function for values around the computed minimum for each estimated parameter in turn. This is helpful to diagnose problems with the optimizer.
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'neighbourhood_size',.5); % width of the window around the computed minimum to be displayed on the diagnostic plots. This width is expressed in percentage deviation. The Inf value is allowed, and will trigger a plot over the entire domain
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'symmetric_plots',true); % ensure that the check plots are symmetric around the minimum. A value of 0 allows to have asymmetric plots, which can be useful if the minimum is close to a domain boundary, or in conjunction with neighbourhood_size = Inf when the domain is not the entire real line
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'number_of_points',20); % number of points around the minimum where the target function is evaluated (for each parameter)
@@ -738,11 +740,6 @@ catch last_error% if check fails, provide info on using calibration if present
rethrow(last_error);
end
-if options_mom_.mode_compute == 0 %We only report value of moments distance at initial value of the parameters
- fprintf('No minimization of moments distance due to ''mode_compute=0''\n')
- return
-end
-
% -------------------------------------------------------------------------
% Step 7a: Method of moments estimation: print some info
% -------------------------------------------------------------------------
@@ -760,27 +757,56 @@ end
if options_mom_.mom.penalized_estimator
fprintf('\n - penalized estimation using deviation from prior mean and weighted with prior precision');
end
-if options_mom_.mode_compute == 1; fprintf('\n - optimizer (mode_compute=1): fmincon');
-elseif options_mom_.mode_compute == 2; fprintf('\n - optimizer (mode_compute=2): continuous simulated annealing');
-elseif options_mom_.mode_compute == 3; fprintf('\n - optimizer (mode_compute=3): fminunc');
-elseif options_mom_.mode_compute == 4; fprintf('\n - optimizer (mode_compute=4): csminwel');
-elseif options_mom_.mode_compute == 5; fprintf('\n - optimizer (mode_compute=5): newrat');
-elseif options_mom_.mode_compute == 6; fprintf('\n - optimizer (mode_compute=6): gmhmaxlik');
-elseif options_mom_.mode_compute == 7; fprintf('\n - optimizer (mode_compute=7): fminsearch');
-elseif options_mom_.mode_compute == 8; fprintf('\n - optimizer (mode_compute=8): Dynare Nelder-Mead simplex');
-elseif options_mom_.mode_compute == 9; fprintf('\n - optimizer (mode_compute=9): CMA-ES');
-elseif options_mom_.mode_compute == 10; fprintf('\n - optimizer (mode_compute=10): simpsa');
-elseif options_mom_.mode_compute == 11; fprintf('\n - optimizer (mode_compute=11): online_auxiliary_filter');
-elseif options_mom_.mode_compute == 12; fprintf('\n - optimizer (mode_compute=12): particleswarm');
-elseif options_mom_.mode_compute == 101; fprintf('\n - optimizer (mode_compute=101): SolveOpt');
-elseif options_mom_.mode_compute == 102; fprintf('\n - optimizer (mode_compute=102): simulannealbnd');
-elseif options_mom_.mode_compute == 13; fprintf('\n - optimizer (mode_compute=13): lsqnonlin');
-elseif ischar(minimizer_algorithm); fprintf(['\n - user-defined optimizer: ' minimizer_algorithm]);
-else
- error('method_of_moments: Unknown optimizer, please contact the developers ')
-end
-if options_mom_.silent_optimizer
- fprintf(' (silent)');
+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';
+ else
+ str = ' (additional_optimizer_steps';
+ end
+ switch optimizer_vec{i}
+ case 0
+ fprintf('\n %s=0): no minimization',str);
+ case 1
+ fprintf('\n %s=1): fmincon',str);
+ case 2
+ fprintf('\n %s=2): continuous simulated annealing',str);
+ case 3
+ fprintf('\n %s=3): fminunc',str);
+ case 4
+ fprintf('\n %s=4): csminwel',str);
+ case 5
+ fprintf('\n %s=5): newrat',str);
+ case 6
+ fprintf('\n %s=6): gmhmaxlik',str);
+ case 7
+ fprintf('\n %s=7): fminsearch',str);
+ case 8
+ fprintf('\n %s=8): Dynare Nelder-Mead simplex',str);
+ case 9
+ fprintf('\n %s=9): CMA-ES',str);
+ case 10
+ fprintf('\n %s=10): simpsa',str);
+ case 11
+ fprintf('\n %s=11): online_auxiliary_filter',str);
+ case 12
+ fprintf('\n %s=12): particleswarm',str);
+ case 101
+ fprintf('\n %s=101): SolveOpt',str);
+ case 102
+ fprintf('\n %s=102): simulannealbnd',str);
+ case 13
+ fprintf('\n %s=13): lsqnonlin',str);
+ otherwise
+ if ischar(optimizer_vec{i})
+ fprintf('\n %s=%s): user-defined',str,optimizer_vec{i});
+ else
+ error('method_of_moments: Unknown optimizer, please contact the developers ')
+ end
+ end
+ if options_mom_.silent_optimizer
+ fprintf(' (silent)');
+ end
end
fprintf('\n - perturbation order: %d', options_mom_.order)
if options_mom_.order > 1 && options_mom_.pruning
@@ -802,8 +828,6 @@ 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
-optimizer_vec=[options_mom_.mode_compute,options_mom_.additional_optimizer_steps]; % at each stage one can possibly use different optimizers sequentially
-
for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
fprintf('Estimation stage %u\n',stage_iter);
Woptflag = false;
@@ -849,15 +873,20 @@ for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
end
for optim_iter= 1:length(optimizer_vec)
- if optimizer_vec(optim_iter)==13
- options_mom_.vector_output = true;
+ 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
- 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, [],...
- Bounds, oo_, estim_params_, M_, options_mom_);
- if options_mom_.vector_output
- fval = fval'*fval;
+ if optimizer_vec{optim_iter}==13
+ options_mom_.vector_output = true;
+ else
+ 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, [],...
+ Bounds, oo_, estim_params_, M_, options_mom_);
+ if options_mom_.vector_output
+ fval = fval'*fval;
+ end
end
fprintf('\nStage %d Iteration %d: value of minimized moment distance objective function: %12.10f.\n',stage_iter,optim_iter,fval)
if options_mom_.mom.verbose
diff --git a/matlab/method_of_moments/method_of_moments_check_plot.m b/matlab/method_of_moments/method_of_moments_check_plot.m
index 24b3ebb61..ceb1a15b7 100644
--- a/matlab/method_of_moments/method_of_moments_check_plot.m
+++ b/matlab/method_of_moments/method_of_moments_check_plot.m
@@ -2,7 +2,7 @@ function method_of_moments_check_plot(fun,xparam,SE_vec,options_,M_,estim_params
% Checks the estimated local minimum of the moment's distance objective
-% Copyright (C) 2020 Dynare Team
+% Copyright (C) 2020-2021 Dynare Team
%
% This file is part of Dynare.
%
diff --git a/matlab/method_of_moments/method_of_moments_data_moments.m b/matlab/method_of_moments/method_of_moments_data_moments.m
index 15e756955..26bf7a4dc 100644
--- a/matlab/method_of_moments/method_of_moments_data_moments.m
+++ b/matlab/method_of_moments/method_of_moments_data_moments.m
@@ -16,7 +16,7 @@ function [dataMoments, m_data] = method_of_moments_data_moments(data, oo_, match
% o method_of_moments.m
% o method_of_moments_objective_function.m
% =========================================================================
-% Copyright (C) 2020 Dynare Team
+% Copyright (C) 2020-2021 Dynare Team
%
% This file is part of Dynare.
%
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 b97c7725f..745f8570e 100644
--- a/matlab/method_of_moments/method_of_moments_objective_function.m
+++ b/matlab/method_of_moments/method_of_moments_objective_function.m
@@ -31,7 +31,7 @@ function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_] = method_o
% o resol
% o set_all_parameters
% =========================================================================
-% Copyright (C) 2020 Dynare Team
+% Copyright (C) 2020-2021 Dynare Team
%
% This file is part of Dynare.
%
@@ -232,7 +232,7 @@ elseif strcmp(options_mom_.mom.mom_method,'SMM')
i_ME = setdiff([1:size(M_.H,1)],find(diag(M_.H) == 0)); % find ME with 0 variance
chol_S = chol(M_.H(i_ME,i_ME)); %decompose rest
shock_mat=zeros(size(options_mom_.mom.ME_shock_series)); %initialize
- shock_mat(:,i_ME)=options_mom_.mom.ME_shock_series(:,i_exo_var)*chol_S;
+ shock_mat(:,i_ME)=options_mom_.mom.ME_shock_series(:,i_ME)*chol_S;
y_sim = y_sim+shock_mat;
end
diff --git a/matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m b/matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m
index 7dde93568..ad2db20f7 100644
--- a/matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m
+++ b/matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m
@@ -19,7 +19,7 @@ function W_opt = method_of_moments_optimal_weighting_matrix(m_data, moments, q_l
% This function calls:
% o CorrMatrix (embedded)
% =========================================================================
-% Copyright (C) 2020 Dynare Team
+% Copyright (C) 2020-2021 Dynare Team
%
% This file is part of Dynare.
%
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 7bfd321b5..459f583f7 100644
--- a/matlab/method_of_moments/method_of_moments_standard_errors.m
+++ b/matlab/method_of_moments/method_of_moments_standard_errors.m
@@ -29,7 +29,7 @@ function [SE_values, Asympt_Var] = method_of_moments_standard_errors(xparam, obj
% o SMM_objective_function.m
% o method_of_moments_optimal_weighting_matrix
% =========================================================================
-% Copyright (C) 2020 Dynare Team
+% Copyright (C) 2020-2021 Dynare Team
%
% This file is part of Dynare.
%
diff --git a/tests/.gitignore b/tests/.gitignore
index 722a27d5d..275d8736e 100644
--- a/tests/.gitignore
+++ b/tests/.gitignore
@@ -50,8 +50,10 @@ wsOct
!/ep/mean_preserving_spread.m
!/ep/rbcii_steady_state.m
!/estimation/fsdat_simul.m
-!/estimation/method_of_moments/RBC_MoM_steady_helper.m
-!/estimation/method_of_moments/RBC_Andreasen_Data_2.mat
+!/estimation/method_of_moments/RBC/RBC_MoM_steady_helper.m
+!/estimation/method_of_moments/RBC/RBC_Andreasen_Data_2.mat
+!/estimation/method_of_moments/AFVRR/AFVRR_data.mat
+!/estimation/method_of_moments/AFVRR/AFVRR_steady_helper.m
!/expectations/expectation_ss_old_steadystate.m
!/external_function/extFunDeriv.m
!/external_function/extFunNoDerivs.m
diff --git a/tests/Makefile.am b/tests/Makefile.am
index fd050916a..cd60e0473 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -50,10 +50,14 @@ MODFILES = \
estimation/MH_recover/fs2000_recover_3.mod \
estimation/t_proposal/fs2000_student.mod \
estimation/tune_mh_jscale/fs2000.mod \
- estimation/method_of_moments/AnScho_MoM.mod \
- estimation/method_of_moments/RBC_MoM_Andreasen.mod \
- estimation/method_of_moments/RBC_MoM_SMM_ME.mod \
- estimation/method_of_moments/RBC_MoM_prefilter.mod \
+ estimation/method_of_moments/AnScho/AnScho_MoM.mod \
+ estimation/method_of_moments/RBC/RBC_MoM_Andreasen.mod \
+ 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/AFVRR/AFVRR_M0.mod \
+ estimation/method_of_moments/AFVRR/AFVRR_MFB.mod \
+ estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod \
moments/example1_var_decomp.mod \
moments/example1_bp_test.mod \
moments/test_AR1_spectral_density.mod \
@@ -835,6 +839,10 @@ particle: m/particle o/particle
m/particle: $(patsubst %.mod, %.m.trs, $(PARTICLEFILES))
o/particle: $(patsubst %.mod, %.o.trs, $(PARTICLEFILES))
+method_of_moments: m/method_of_moments o/method_of_moments
+m/method_of_moments: $(patsubst %.mod, %.m.trs, $(filter estimation/method_of_moments/%.mod, $(MODFILES)))
+o/method_of_moments: $(patsubst %.mod, %.o.trs, $(filter estimation/method_of_moments/%.mod, $(MODFILES)))
+
# Matlab TRS Files
M_TRS_FILES = $(patsubst %.mod, %.m.trs, $(MODFILES))
M_TRS_FILES += run_block_byte_tests_matlab.m.trs \
@@ -984,8 +992,10 @@ EXTRA_DIST = \
lmmcp/sw-common-header.inc \
lmmcp/sw-common-footer.inc \
estimation/tune_mh_jscale/fs2000.inc \
- estimation/method_of_moments/RBC_MoM_common.inc \
- estimation/method_of_moments/RBC_MoM_steady_helper.m \
+ estimation/method_of_moments/RBC/RBC_MoM_common.inc \
+ estimation/method_of_moments/RBC/RBC_MoM_steady_helper.m \
+ estimation/method_of_moments/AFVRR/AFVRR_common.inc \
+ estimation/method_of_moments/AFVRR/AFVRR_steady_helper.m \
histval_initval_file_unit_tests.m \
histval_initval_file/my_assert.m \
histval_initval_file/ramst_data.xls \
diff --git a/tests/estimation/method_of_moments/AFVRR/AFVRR_M0.mod b/tests/estimation/method_of_moments/AFVRR/AFVRR_M0.mod
new file mode 100644
index 000000000..8e51ac513
--- /dev/null
+++ b/tests/estimation/method_of_moments/AFVRR/AFVRR_M0.mod
@@ -0,0 +1,299 @@
+% DSGE model based on replication files of
+% Andreasen, Fernandez-Villaverde, Rubio-Ramirez (2018), The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications, Review of Economic Studies, 85, p. 1-49
+% Adapted for Dynare by Willi Mutschler (@wmutschl, willi@mutschler.eu), Jan 2021
+% =========================================================================
+% 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 .
+% =========================================================================
+
+% This is the benchmark model with no feedback M_0
+% Original code RunGMM_standardModel_RRA.m by Martin M. Andreasen, Jan 2016
+
+@#include "AFVRR_common.inc"
+
+%--------------------------------------------------------------------------
+% Parameter calibration taken from RunGMM_standardModel_RRA.m
+%--------------------------------------------------------------------------
+% fixed parameters
+INHABIT = 1;
+PHI1 = 4;
+PHI4 = 1;
+KAPAone = 0;
+DELTA = 0.025;
+THETA = 0.36;
+ETA = 6;
+CHI = 0;
+CONSxhr40 = 0;
+BETTAxhr = 0;
+BETTAxhr40= 0;
+RHOD = 0;
+GAMA = 0.9999;
+CONSxhr20 = 0;
+
+% estimated parameters
+BETTA = 0.999544966118000;
+B = 0.668859504661000;
+H = 0.342483445196000;
+PHI2 = 0.997924964981000;
+RRA = 662.7953149595370;
+KAPAtwo = 5.516226495551000;
+ALFA = 0.809462321180000;
+RHOR = 0.643873352513000;
+BETTAPAI = 1.270087844103000;
+BETTAY = 0.031812764291000;
+MYYPS = 1.001189151180000;
+MYZ = 1.005286347928000;
+RHOA = 0.743239127127000;
+RHOG = 0.793929380230000;
+PAI = 1.012163659169000;
+GoY = 0.206594858866000;
+STDA = 0.016586292524000;
+STDG = 0.041220613851000;
+STDD = 0.013534473123000;
+
+% endogenous parameters set via steady state, no need to initialize
+%PHIzero = ;
+%AA = ;
+%PHI3 = ;
+%negVf = ;
+
+model_diagnostics;
+% Model diagnostics show that some parameters are endogenously determined
+% via the steady state, so we run steady to calibrate all parameters
+steady;
+model_diagnostics;
+% Now all parameters are determined
+
+resid;
+check;
+
+%--------------------------------------------------------------------------
+% Shock distribution
+%--------------------------------------------------------------------------
+shocks;
+var eps_a = STDA^2;
+var eps_d = STDD^2;
+var eps_g = STDG^2;
+end;
+
+%--------------------------------------------------------------------------
+% Estimated Params block - these parameters will be estimated, we
+% initialize at calibrated values
+%--------------------------------------------------------------------------
+estimated_params;
+BETTA;
+B;
+H;
+PHI2;
+RRA;
+KAPAtwo;
+ALFA;
+RHOR;
+BETTAPAI;
+BETTAY;
+MYYPS;
+MYZ;
+RHOA;
+RHOG;
+PAI;
+GoY;
+stderr eps_a;
+stderr eps_g;
+stderr eps_d;
+end;
+
+estimated_params_init(use_calibration);
+end;
+
+%--------------------------------------------------------------------------
+% Compare whether toolbox yields equivalent moments at second order
+%--------------------------------------------------------------------------
+% Note that we compare results for orderApp=1|2 and not for orderApp=3, because
+% there is a small error in the replication files of the original article in the
+% computation of the covariance matrix of the extended innovations vector.
+% The authors have been contacted, fixed it, and report that the results
+% change only slightly at orderApp=3 to what they report in the paper. At
+% orderApp=2 all is correct and so the following part tests whether we get
+% the same model moments at the calibrated parameters (we do not optimize).
+% We compare it to the replication file RunGMM_standardModel_RRA.m with the
+% following settings: orderApp=1|2, seOn=0, q_lag=10, weighting=1;
+% scaled=0; optimizer=0; estimator=1; momentSet=2;
+%
+% Output of the replication files for orderApp=1
+AndreasenEtAl.Q1 = 23893.072;
+AndreasenEtAl.moments1 =[ % note that we reshuffeled to be compatible with our matched moments block
+ {[ 1]} {'Ex' } {'Gr_C '} {' ' } {'0.024388' } {'0.023764' }
+ {[ 2]} {'Ex' } {'Gr_I '} {' ' } {'0.031046' } {'0.028517' }
+ {[ 3]} {'Ex' } {'Infl ' } {' ' } {'0.03757' } {'0.048361' }
+ {[ 4]} {'Ex' } {'r1 ' } {' ' } {'0.056048' } {'0.073945' }
+ {[ 5]} {'Ex' } {'r40 ' } {' ' } {'0.069929' } {'0.073945' }
+ {[ 6]} {'Ex' } {'xhr40 '} {' ' } {'0.017237' } {'0' }
+ {[ 7]} {'Ex' } {'GoY '} {' ' } {'-1.5745' } {'-1.577' }
+ {[ 8]} {'Ex' } {'hours '} {' ' } {'-0.043353' } {'-0.042861' }
+ {[ 9]} {'Exx' } {'Gr_C '} {'Gr_C '} {'0.0013159' } {'0.0011816' }
+ {[17]} {'Exx' } {'Gr_C '} {'Gr_I '} {'0.0021789' } {'0.0016052' }
+ {[18]} {'Exx' } {'Gr_C '} {'Infl ' } {'0.00067495' } {'0.00090947' }
+ {[19]} {'Exx' } {'Gr_C '} {'r1 ' } {'0.0011655' } {'0.0016016' }
+ {[20]} {'Exx' } {'Gr_C '} {'r40 ' } {'0.0015906' } {'0.0017076' }
+ {[21]} {'Exx' } {'Gr_C '} {'xhr40 '} {'0.0020911' } {'0.0013997' }
+ {[10]} {'Exx' } {'Gr_I '} {'Gr_I '} {'0.0089104' } {'0.0055317' }
+ {[22]} {'Exx' } {'Gr_I '} {'Infl ' } {'0.00063139' } {'0.00050106' }
+ {[23]} {'Exx' } {'Gr_I '} {'r1 ' } {'0.0011031' } {'0.0018178' }
+ {[24]} {'Exx' } {'Gr_I '} {'r40 ' } {'0.0018445' } {'0.0020186' }
+ {[25]} {'Exx' } {'Gr_I '} {'xhr40 '} {'0.00095556' } {'0.0064471' }
+ {[11]} {'Exx' } {'Infl ' } {'Infl ' } {'0.0020268' } {'0.0030519' }
+ {[26]} {'Exx' } {'Infl ' } {'r1 ' } {'0.0025263' } {'0.0042181' }
+ {[27]} {'Exx' } {'Infl ' } {'r40 ' } {'0.0029126' } {'0.0039217' }
+ {[28]} {'Exx' } {'Infl ' } {'xhr40 '} {'-0.00077101'} {'-0.0019975' }
+ {[12]} {'Exx' } {'r1 ' } {'r1 ' } {'0.0038708' } {'0.0061403' }
+ {[29]} {'Exx' } {'r1 ' } {'r40 ' } {'0.0044773' } {'0.0058343' }
+ {[30]} {'Exx' } {'r1 ' } {'xhr40 '} {'-0.00048202'} {'-0.00089501'}
+ {[13]} {'Exx' } {'r40 ' } {'r40 ' } {'0.0054664' } {'0.0056883' }
+ {[31]} {'Exx' } {'r40 ' } {'xhr40 '} {'0.00053864' } {'-0.00041184'}
+ {[14]} {'Exx' } {'xhr40 '} {'xhr40 '} {'0.053097' } {'0.016255' }
+ {[15]} {'Exx' } {'GoY '} {'GoY '} {'2.4863' } {'2.4919' }
+ {[16]} {'Exx' } {'hours '} {'hours '} {'0.0018799' } {'0.0018384' }
+ {[32]} {'Exx1'} {'Gr_C '} {'Gr_C '} {'0.00077917' } {'0.00065543' }
+ {[33]} {'Exx1'} {'Gr_I '} {'Gr_I '} {'0.0050104' } {'0.0033626' }
+ {[34]} {'Exx1'} {'Infl ' } {'Infl ' } {'0.0019503' } {'0.0029033' }
+ {[35]} {'Exx1'} {'r1 ' } {'r1 ' } {'0.0038509' } {'0.006112' }
+ {[36]} {'Exx1'} {'r40 ' } {'r40 ' } {'0.0054699' } {'0.005683' }
+ {[37]} {'Exx1'} {'xhr40 '} {'xhr40 '} {'-0.00098295'} {'3.3307e-16' }
+ {[38]} {'Exx1'} {'GoY '} {'GoY '} {'2.4868' } {'2.4912' }
+ {[39]} {'Exx1'} {'hours '} {'hours '} {'0.0018799' } {'0.0018378' }
+];
+
+% Output of the replication files for orderApp=2
+AndreasenEtAl.Q2 = 65.8269;
+AndreasenEtAl.moments2 =[ % note that we reshuffeled to be compatible with our matched moments block
+ {[ 1]} {'Ex' } {'Gr_C '} {' ' } {'0.024388' } {'0.023764' }
+ {[ 2]} {'Ex' } {'Gr_I '} {' ' } {'0.031046' } {'0.028517' }
+ {[ 3]} {'Ex' } {'Infl ' } {' ' } {'0.03757' } {'0.034882' }
+ {[ 4]} {'Ex' } {'r1 ' } {' ' } {'0.056048' } {'0.056542' }
+ {[ 5]} {'Ex' } {'r40 ' } {' ' } {'0.069929' } {'0.070145' }
+ {[ 6]} {'Ex' } {'xhr40 '} {' ' } {'0.017237' } {'0.020825' }
+ {[ 7]} {'Ex' } {'GoY '} {' ' } {'-1.5745' } {'-1.5748' }
+ {[ 8]} {'Ex' } {'hours '} {' ' } {'-0.043353' } {'-0.04335' }
+ {[ 9]} {'Exx' } {'Gr_C '} {'Gr_C '} {'0.0013159' } {'0.001205' }
+ {[17]} {'Exx' } {'Gr_C '} {'Gr_I '} {'0.0021789' } {'0.0016067' }
+ {[18]} {'Exx' } {'Gr_C '} {'Infl ' } {'0.00067495' } {'0.00059406'}
+ {[19]} {'Exx' } {'Gr_C '} {'r1 ' } {'0.0011655' } {'0.0011949' }
+ {[20]} {'Exx' } {'Gr_C '} {'r40 ' } {'0.0015906' } {'0.0016104' }
+ {[21]} {'Exx' } {'Gr_C '} {'xhr40 '} {'0.0020911' } {'0.0020245' }
+ {[10]} {'Exx' } {'Gr_I '} {'Gr_I '} {'0.0089104' } {'0.0060254' }
+ {[22]} {'Exx' } {'Gr_I '} {'Infl ' } {'0.00063139' } {'8.3563e-05'}
+ {[23]} {'Exx' } {'Gr_I '} {'r1 ' } {'0.0011031' } {'0.0013176' }
+ {[24]} {'Exx' } {'Gr_I '} {'r40 ' } {'0.0018445' } {'0.0019042' }
+ {[25]} {'Exx' } {'Gr_I '} {'xhr40 '} {'0.00095556' } {'0.0064261' }
+ {[11]} {'Exx' } {'Infl ' } {'Infl ' } {'0.0020268' } {'0.0020735' }
+ {[26]} {'Exx' } {'Infl ' } {'r1 ' } {'0.0025263' } {'0.0027621' }
+ {[27]} {'Exx' } {'Infl ' } {'r40 ' } {'0.0029126' } {'0.0029257' }
+ {[28]} {'Exx' } {'Infl ' } {'xhr40 '} {'-0.00077101'} {'-0.0012165'}
+ {[12]} {'Exx' } {'r1 ' } {'r1 ' } {'0.0038708' } {'0.0040235' }
+ {[29]} {'Exx' } {'r1 ' } {'r40 ' } {'0.0044773' } {'0.0044702' }
+ {[30]} {'Exx' } {'r1 ' } {'xhr40 '} {'-0.00048202'} {'0.00030542'}
+ {[13]} {'Exx' } {'r40 ' } {'r40 ' } {'0.0054664' } {'0.0052718' }
+ {[31]} {'Exx' } {'r40 ' } {'xhr40 '} {'0.00053864' } {'0.0010045' }
+ {[14]} {'Exx' } {'xhr40 '} {'xhr40 '} {'0.053097' } {'0.018416' }
+ {[15]} {'Exx' } {'GoY '} {'GoY '} {'2.4863' } {'2.4853' }
+ {[16]} {'Exx' } {'hours '} {'hours '} {'0.0018799' } {'0.0018806' }
+ {[32]} {'Exx1'} {'Gr_C '} {'Gr_C '} {'0.00077917' } {'0.00067309'}
+ {[33]} {'Exx1'} {'Gr_I '} {'Gr_I '} {'0.0050104' } {'0.0033293' }
+ {[34]} {'Exx1'} {'Infl ' } {'Infl ' } {'0.0019503' } {'0.0019223' }
+ {[35]} {'Exx1'} {'r1 ' } {'r1 ' } {'0.0038509' } {'0.0039949' }
+ {[36]} {'Exx1'} {'r40 ' } {'r40 ' } {'0.0054699' } {'0.0052659' }
+ {[37]} {'Exx1'} {'xhr40 '} {'xhr40 '} {'-0.00098295'} {'0.0004337' }
+ {[38]} {'Exx1'} {'GoY '} {'GoY '} {'2.4868' } {'2.4846' }
+ {[39]} {'Exx1'} {'hours '} {'hours '} {'0.0018799' } {'0.00188' }
+];
+
+@#for orderApp in 1:2
+
+method_of_moments(
+ mom_method = GMM % method of moments method; possible values: GMM|SMM
+ , datafile = 'AFVRR_data.mat' % name of filename with data
+ , bartlett_kernel_lag = 10 % bandwith in optimal weighting matrix
+ , order = @{orderApp} % order of Taylor approximation in perturbation
+ , pruning % use pruned state space system at higher-order
+ % , verbose % display and store intermediate estimation results
+ , weighting_matrix = ['DIAGONAL'] % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
+ % , TeX % print TeX tables and graphics
+ % Optimization options that can be set by the user in the mod file, otherwise default values are provided
+ %, huge_number=1D10 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+ , mode_compute = 0 % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
+ , optim = ('TolFun', 1e-6
+ ,'TolX', 1e-6
+ ,'MaxIter', 3000
+ ,'MaxFunEvals', 1D6
+ ,'UseParallel' , 1
+ %,'Jacobian' , '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_standard_errors
+ , se_tolx=1e-10
+);
+
+% Check results
+
+fprintf('****************************************************************\n')
+fprintf('Compare Results for perturbation order @{orderApp}\n')
+fprintf('****************************************************************\n')
+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;
+
+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))])
+
+if norm(dev_modelmoments)> 1e-4
+ error('Something wrong in the computation of moments at order @{orderApp}')
+end
+
+@#endfor
+
+%--------------------------------------------------------------------------
+% Replicate estimation at orderApp=3
+%--------------------------------------------------------------------------
+@#ifdef DoEstimation
+method_of_moments(
+ mom_method = GMM % method of moments method; possible values: GMM|SMM
+ , datafile = 'AFVRR_data.mat' % name of filename with data
+ , bartlett_kernel_lag = 10 % bandwith in optimal weighting matrix
+ , order = 3 % order of Taylor approximation in perturbation
+ , pruning % use pruned state space system at higher-order
+ % , verbose % display and store intermediate estimation results
+ , weighting_matrix = ['DIAGONAL', 'OPTIMAL'] % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
+ % , TeX % print TeX tables and graphics
+ % Optimization options that can be set by the user in the mod file, otherwise default values are provided
+ %, huge_number=1D10 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+ , mode_compute = 13 % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
+ , additional_optimizer_steps = [13]
+ , optim = ('TolFun', 1e-6
+ ,'TolX', 1e-6
+ ,'MaxIter', 3000
+ ,'MaxFunEvals', 1D6
+ ,'UseParallel' , 1
+ %,'Jacobian' , '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_standard_errors
+ , se_tolx=1e-10
+);
+@#endif
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/AFVRR/AFVRR_MFB.mod b/tests/estimation/method_of_moments/AFVRR/AFVRR_MFB.mod
new file mode 100644
index 000000000..450739ad3
--- /dev/null
+++ b/tests/estimation/method_of_moments/AFVRR/AFVRR_MFB.mod
@@ -0,0 +1,300 @@
+% DSGE model based on replication files of
+% Andreasen, Fernandez-Villaverde, Rubio-Ramirez (2018), The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications, Review of Economic Studies, 85, p. 1-49
+% Adapted for Dynare by Willi Mutschler (@wmutschl, willi@mutschler.eu), Jan 2021
+% =========================================================================
+% 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 .
+% =========================================================================
+
+% This is the model with Feedback M_FB
+% Original code RunGMM_Feedback_estim_RRA.m by Martin M. Andreasen, Jan 2016
+
+@#include "AFVRR_common.inc"
+
+%--------------------------------------------------------------------------
+% Parameter calibration taken from RunGMM_Feedback_estim_RRA.m
+%--------------------------------------------------------------------------
+% fixed parameters
+INHABIT = 1;
+PHI1 = 4;
+PHI4 = 1;
+KAPAone = 0;
+DELTA = 0.025;
+THETA = 0.36;
+ETA = 6;
+CHI = 0;
+BETTAxhr = 0;
+BETTAxhr40= 0;
+RHOD = 0;
+GAMA = 0.9999;
+CONSxhr20 = 0;
+
+% estimated parameters
+BETTA = 0.997007023687000;
+B = 0.692501768577000;
+H = 0.339214495653000;
+PHI2 = 0.688555040951000;
+RRA = 24.346514272871001;
+KAPAtwo = 10.018421876923000;
+ALFA = 0.792507553312000;
+RHOR = 0.849194030384000;
+BETTAPAI = 2.060579322980000;
+BETTAY = 0.220573712342000;
+MYYPS = 1.001016690133000;
+MYZ = 1.005356313981000;
+RHOA = 0.784141391843000;
+RHOG = 0.816924540497000;
+PAI = 1.011924196487000;
+CONSxhr40 = 0.878774662208000;
+GoY = 0.207110300602000;
+STDA = 0.013024450606000;
+STDG = 0.051049871928000;
+STDD = 0.008877423780000;
+
+% endogenous parameters set via steady state, no need to initialize
+%PHIzero = ;
+%AA = ;
+%PHI3 = ;
+%negVf = ;
+
+model_diagnostics;
+% Model diagnostics show that some parameters are endogenously determined
+% via the steady state, so we run steady to calibrate all parameters
+steady;
+model_diagnostics;
+% Now all parameters are determined
+
+resid;
+check;
+
+%--------------------------------------------------------------------------
+% Shock distribution
+%--------------------------------------------------------------------------
+shocks;
+var eps_a = STDA^2;
+var eps_d = STDD^2;
+var eps_g = STDG^2;
+end;
+
+%--------------------------------------------------------------------------
+% Estimated Params block - these parameters will be estimated, we
+% initialize at calibrated values
+%--------------------------------------------------------------------------
+estimated_params;
+BETTA;
+B;
+H;
+PHI2;
+RRA;
+KAPAtwo;
+ALFA;
+RHOR;
+BETTAPAI;
+BETTAY;
+MYYPS;
+MYZ;
+RHOA;
+RHOG;
+PAI;
+CONSxhr40;
+GoY;
+stderr eps_a;
+stderr eps_g;
+stderr eps_d;
+end;
+
+estimated_params_init(use_calibration);
+end;
+
+%--------------------------------------------------------------------------
+% Compare whether toolbox yields equivalent moments at second order
+%--------------------------------------------------------------------------
+% Note that we compare results for orderApp=1|2 and not for orderApp=3, because
+% there is a small error in the replication files of the original article in the
+% computation of the covariance matrix of the extended innovations vector.
+% The authors have been contacted, fixed it, and report that the results
+% change only slightly at orderApp=3 to what they report in the paper. At
+% orderApp=2 all is correct and so the following part tests whether we get
+% the same model moments at the calibrated parameters (we do not optimize).
+% We compare it to the replication file RunGMM_Feedback_estim_RRA.m with the
+% following settings: orderApp=1|2, seOn=0, q_lag=10, weighting=1;
+% scaled=0; optimizer=0; estimator=1; momentSet=2;
+%
+% Output of the replication files for orderApp=1
+AndreasenEtAl.Q1 = 201778.9697;
+AndreasenEtAl.moments1 =[ % note that we reshuffeled to be compatible with our matched moments block
+ {[ 1]} {'Ex' } {'Gr_C '} {' ' } {'0.024388' } {'0.023654' }
+ {[ 2]} {'Ex' } {'Gr_I '} {' ' } {'0.031046' } {'0.027719' }
+ {[ 3]} {'Ex' } {'Infl ' } {' ' } {'0.03757' } {'0.047415' }
+ {[ 4]} {'Ex' } {'r1 ' } {' ' } {'0.056048' } {'0.083059' }
+ {[ 5]} {'Ex' } {'r40 ' } {' ' } {'0.069929' } {'0.083059' }
+ {[ 6]} {'Ex' } {'xhr40 '} {' ' } {'0.017237' } {'0' }
+ {[ 7]} {'Ex' } {'GoY '} {' ' } {'-1.5745' } {'-1.5745' }
+ {[ 8]} {'Ex' } {'hours '} {' ' } {'-0.043353' } {'-0.043245' }
+ {[ 9]} {'Exx' } {'Gr_C '} {'Gr_C '} {'0.0013159' } {'0.0012253' }
+ {[17]} {'Exx' } {'Gr_C '} {'Gr_I '} {'0.0021789' } {'0.0015117' }
+ {[18]} {'Exx' } {'Gr_C '} {'Infl ' } {'0.00067495' } {'0.00080078' }
+ {[19]} {'Exx' } {'Gr_C '} {'r1 ' } {'0.0011655' } {'0.00182' }
+ {[20]} {'Exx' } {'Gr_C '} {'r40 ' } {'0.0015906' } {'0.001913' }
+ {[21]} {'Exx' } {'Gr_C '} {'xhr40 '} {'0.0020911' } {'0.0016326' }
+ {[10]} {'Exx' } {'Gr_I '} {'Gr_I '} {'0.0089104' } {'0.0040112' }
+ {[22]} {'Exx' } {'Gr_I '} {'Infl ' } {'0.00063139' } {'0.00060604' }
+ {[23]} {'Exx' } {'Gr_I '} {'r1 ' } {'0.0011031' } {'0.0021426' }
+ {[24]} {'Exx' } {'Gr_I '} {'r40 ' } {'0.0018445' } {'0.0022348' }
+ {[25]} {'Exx' } {'Gr_I '} {'xhr40 '} {'0.00095556' } {'0.0039852' }
+ {[11]} {'Exx' } {'Infl ' } {'Infl ' } {'0.0020268' } {'0.0030058' }
+ {[26]} {'Exx' } {'Infl ' } {'r1 ' } {'0.0025263' } {'0.0044951' }
+ {[27]} {'Exx' } {'Infl ' } {'r40 ' } {'0.0029126' } {'0.0042225' }
+ {[28]} {'Exx' } {'Infl ' } {'xhr40 '} {'-0.00077101'} {'-0.0021222' }
+ {[12]} {'Exx' } {'r1 ' } {'r1 ' } {'0.0038708' } {'0.0074776' }
+ {[29]} {'Exx' } {'r1 ' } {'r40 ' } {'0.0044773' } {'0.0071906' }
+ {[30]} {'Exx' } {'r1 ' } {'xhr40 '} {'-0.00048202'} {'-0.0006736' }
+ {[13]} {'Exx' } {'r40 ' } {'r40 ' } {'0.0054664' } {'0.0070599' }
+ {[31]} {'Exx' } {'r40 ' } {'xhr40 '} {'0.00053864' } {'-0.00036735'}
+ {[14]} {'Exx' } {'xhr40 '} {'xhr40 '} {'0.053097' } {'0.014516' }
+ {[15]} {'Exx' } {'GoY '} {'GoY '} {'2.4863' } {'2.4866' }
+ {[16]} {'Exx' } {'hours '} {'hours '} {'0.0018799' } {'0.0018713' }
+ {[32]} {'Exx1'} {'Gr_C '} {'Gr_C '} {'0.00077917' } {'0.00076856' }
+ {[33]} {'Exx1'} {'Gr_I '} {'Gr_I '} {'0.0050104' } {'0.002163' }
+ {[34]} {'Exx1'} {'Infl ' } {'Infl ' } {'0.0019503' } {'0.0028078' }
+ {[35]} {'Exx1'} {'r1 ' } {'r1 ' } {'0.0038509' } {'0.0074583' }
+ {[36]} {'Exx1'} {'r40 ' } {'r40 ' } {'0.0054699' } {'0.0070551' }
+ {[37]} {'Exx1'} {'xhr40 '} {'xhr40 '} {'-0.00098295'} {'7.2164e-16' }
+ {[38]} {'Exx1'} {'GoY '} {'GoY '} {'2.4868' } {'2.4856' }
+ {[39]} {'Exx1'} {'hours '} {'hours '} {'0.0018799' } {'0.0018708' }
+];
+
+% Output of the replication files for orderApp=2
+AndreasenEtAl.Q2 = 59.3323;
+AndreasenEtAl.moments2 =[ % note that we reshuffeled to be compatible with our matched moments block
+ {[ 1]} {'Ex' } {'Gr_C '} {' ' } {'0.024388' } {'0.023654' }
+ {[ 2]} {'Ex' } {'Gr_I '} {' ' } {'0.031046' } {'0.027719' }
+ {[ 3]} {'Ex' } {'Infl ' } {' ' } {'0.03757' } {'0.034565' }
+ {[ 4]} {'Ex' } {'r1 ' } {' ' } {'0.056048' } {'0.056419' }
+ {[ 5]} {'Ex' } {'r40 ' } {' ' } {'0.069929' } {'0.07087' }
+ {[ 6]} {'Ex' } {'xhr40 '} {' ' } {'0.017237' } {'0.01517' }
+ {[ 7]} {'Ex' } {'GoY '} {' ' } {'-1.5745' } {'-1.5743' }
+ {[ 8]} {'Ex' } {'hours '} {' ' } {'-0.043353' } {'-0.043352' }
+ {[ 9]} {'Exx' } {'Gr_C '} {'Gr_C '} {'0.0013159' } {'0.0012464' }
+ {[17]} {'Exx' } {'Gr_C '} {'Gr_I '} {'0.0021789' } {'0.0015247' }
+ {[18]} {'Exx' } {'Gr_C '} {'Infl ' } {'0.00067495' } {'0.0004867' }
+ {[19]} {'Exx' } {'Gr_C '} {'r1 ' } {'0.0011655' } {'0.0011867' }
+ {[20]} {'Exx' } {'Gr_C '} {'r40 ' } {'0.0015906' } {'0.0016146' }
+ {[21]} {'Exx' } {'Gr_C '} {'xhr40 '} {'0.0020911' } {'0.0021395' }
+ {[10]} {'Exx' } {'Gr_I '} {'Gr_I '} {'0.0089104' } {'0.0043272' }
+ {[22]} {'Exx' } {'Gr_I '} {'Infl ' } {'0.00063139' } {'0.00021752'}
+ {[23]} {'Exx' } {'Gr_I '} {'r1 ' } {'0.0011031' } {'0.0013919' }
+ {[24]} {'Exx' } {'Gr_I '} {'r40 ' } {'0.0018445' } {'0.0018899' }
+ {[25]} {'Exx' } {'Gr_I '} {'xhr40 '} {'0.00095556' } {'0.0037854' }
+ {[11]} {'Exx' } {'Infl ' } {'Infl ' } {'0.0020268' } {'0.0021043' }
+ {[26]} {'Exx' } {'Infl ' } {'r1 ' } {'0.0025263' } {'0.0026571' }
+ {[27]} {'Exx' } {'Infl ' } {'r40 ' } {'0.0029126' } {'0.0028566' }
+ {[28]} {'Exx' } {'Infl ' } {'xhr40 '} {'-0.00077101'} {'-0.0016279'}
+ {[12]} {'Exx' } {'r1 ' } {'r1 ' } {'0.0038708' } {'0.0039136' }
+ {[29]} {'Exx' } {'r1 ' } {'r40 ' } {'0.0044773' } {'0.0044118' }
+ {[30]} {'Exx' } {'r1 ' } {'xhr40 '} {'-0.00048202'} {'0.00016791'}
+ {[13]} {'Exx' } {'r40 ' } {'r40 ' } {'0.0054664' } {'0.0052851' }
+ {[31]} {'Exx' } {'r40 ' } {'xhr40 '} {'0.00053864' } {'0.00062143'}
+ {[14]} {'Exx' } {'xhr40 '} {'xhr40 '} {'0.053097' } {'0.018126' }
+ {[15]} {'Exx' } {'GoY '} {'GoY '} {'2.4863' } {'2.4863' }
+ {[16]} {'Exx' } {'hours '} {'hours '} {'0.0018799' } {'0.0018806' }
+ {[32]} {'Exx1'} {'Gr_C '} {'Gr_C '} {'0.00077917' } {'0.00078586'}
+ {[33]} {'Exx1'} {'Gr_I '} {'Gr_I '} {'0.0050104' } {'0.0021519' }
+ {[34]} {'Exx1'} {'Infl ' } {'Infl ' } {'0.0019503' } {'0.0019046' }
+ {[35]} {'Exx1'} {'r1 ' } {'r1 ' } {'0.0038509' } {'0.0038939' }
+ {[36]} {'Exx1'} {'r40 ' } {'r40 ' } {'0.0054699' } {'0.0052792' }
+ {[37]} {'Exx1'} {'xhr40 '} {'xhr40 '} {'-0.00098295'} {'0.00023012'}
+ {[38]} {'Exx1'} {'GoY '} {'GoY '} {'2.4868' } {'2.4852' }
+ {[39]} {'Exx1'} {'hours '} {'hours '} {'0.0018799' } {'0.0018801' }
+];
+
+@#for orderApp in 1:2
+
+method_of_moments(
+ mom_method = GMM % method of moments method; possible values: GMM|SMM
+ , datafile = 'AFVRR_data.mat' % name of filename with data
+ , bartlett_kernel_lag = 10 % bandwith in optimal weighting matrix
+ , order = @{orderApp} % order of Taylor approximation in perturbation
+ , pruning % use pruned state space system at higher-order
+ % , verbose % display and store intermediate estimation results
+ , weighting_matrix = ['DIAGONAL'] % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
+ % , TeX % print TeX tables and graphics
+ % Optimization options that can be set by the user in the mod file, otherwise default values are provided
+ %, huge_number=1D10 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+ , mode_compute = 0 % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
+ , optim = ('TolFun', 1e-6
+ ,'TolX', 1e-6
+ ,'MaxIter', 3000
+ ,'MaxFunEvals', 1D6
+ ,'UseParallel' , 1
+ %,'Jacobian' , '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_standard_errors
+ , se_tolx=1e-10
+);
+
+% Check results
+
+fprintf('****************************************************************\n')
+fprintf('Compare Results for perturbation order @{orderApp}\n')
+fprintf('****************************************************************\n')
+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;
+
+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))])
+
+if norm(dev_modelmoments)> 1e-4
+ warning('Something wrong in the computation of moments at order @{orderApp}')
+end
+
+@#endfor
+
+%--------------------------------------------------------------------------
+% Replicate estimation at orderApp=3
+%--------------------------------------------------------------------------
+@#ifdef DoEstimation
+method_of_moments(
+ mom_method = GMM % method of moments method; possible values: GMM|SMM
+ , datafile = 'AFVRR_data.mat' % name of filename with data
+ , bartlett_kernel_lag = 10 % bandwith in optimal weighting matrix
+ , order = 3 % order of Taylor approximation in perturbation
+ , pruning % use pruned state space system at higher-order
+ % , verbose % display and store intermediate estimation results
+ , weighting_matrix = ['DIAGONAL', 'Optimal'] % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
+ % , TeX % print TeX tables and graphics
+ % Optimization options that can be set by the user in the mod file, otherwise default values are provided
+ %, huge_number=1D10 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+ , mode_compute = 13 % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
+ , additional_optimizer_steps = [13]
+ , optim = ('TolFun', 1e-6
+ ,'TolX', 1e-6
+ ,'MaxIter', 3000
+ ,'MaxFunEvals', 1D6
+ ,'UseParallel' , 1
+ %,'Jacobian' , '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_standard_errors
+ , se_tolx=1e-10
+);
+@#endif
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod b/tests/estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod
new file mode 100644
index 000000000..9c069d3a3
--- /dev/null
+++ b/tests/estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod
@@ -0,0 +1,299 @@
+% DSGE model based on replication files of
+% Andreasen, Fernandez-Villaverde, Rubio-Ramirez (2018), The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications, Review of Economic Studies, 85, p. 1-49
+% Adapted for Dynare by Willi Mutschler (@wmutschl, willi@mutschler.eu), Jan 2021
+% =========================================================================
+% 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 .
+% =========================================================================
+
+% This is the model with feedback and calibrated RRA
+% Original code RunGMM_Feedback_estim_RRA_5.m by Martin M. Andreasen, Jan 2016
+
+@#include "AFVRR_common.inc"
+
+%--------------------------------------------------------------------------
+% Parameter calibration taken from RunGMM_Feedback_estim_RRA_5.m
+%--------------------------------------------------------------------------
+% fixed parameters
+INHABIT = 1;
+PHI1 = 4;
+PHI4 = 1;
+KAPAone = 0;
+DELTA = 0.025;
+THETA = 0.36;
+ETA = 6;
+CHI = 0;
+BETTAxhr = 0;
+BETTAxhr40= 0;
+RHOD = 0;
+GAMA = 0.9999;
+CONSxhr20 = 0;
+RRA = 5;
+
+% estimated parameters
+BETTA = 0.996850651147000;
+B = 0.684201133923000;
+H = 0.338754441432000;
+PHI2 = 0.738293581320000;
+KAPAtwo = 11.664785970704999;
+ALFA = 0.831836572237000;
+RHOR = 0.772754520116000;
+BETTAPAI = 3.020381242896000;
+BETTAY = 0.288367683973000;
+MYYPS = 1.000911709188000;
+MYZ = 1.005433723022000;
+RHOA = 0.749465413198000;
+RHOG = 0.847225569814000;
+PAI = 1.010428794858000;
+CONSxhr40 = 0.992863217133000;
+GoY = 0.207099399789000;
+STDA = 0.015621059978000;
+STDG = 0.047539390956000;
+STDD = 0.008623441943000;
+
+% endogenous parameters set via steady state, no need to initialize
+%PHIzero = ;
+%AA = ;
+%PHI3 = ;
+%negVf = ;
+
+model_diagnostics;
+% Model diagnostics show that some parameters are endogenously determined
+% via the steady state, so we run steady to calibrate all parameters
+steady;
+model_diagnostics;
+% Now all parameters are determined
+
+resid;
+check;
+
+%--------------------------------------------------------------------------
+% Shock distribution
+%--------------------------------------------------------------------------
+shocks;
+var eps_a = STDA^2;
+var eps_d = STDD^2;
+var eps_g = STDG^2;
+end;
+
+%--------------------------------------------------------------------------
+% Estimated Params block - these parameters will be estimated, we
+% initialize at calibrated values
+%--------------------------------------------------------------------------
+estimated_params;
+BETTA;
+B;
+H;
+PHI2;
+KAPAtwo;
+ALFA;
+RHOR;
+BETTAPAI;
+BETTAY;
+MYYPS;
+MYZ;
+RHOA;
+RHOG;
+PAI;
+CONSxhr40;
+GoY;
+stderr eps_a;
+stderr eps_g;
+stderr eps_d;
+end;
+
+estimated_params_init(use_calibration);
+end;
+
+%--------------------------------------------------------------------------
+% Compare whether toolbox yields equivalent moments at second order
+%--------------------------------------------------------------------------
+% Note that we compare results for orderApp=1|2 and not for orderApp=3, because
+% there is a small error in the replication files of the original article in the
+% computation of the covariance matrix of the extended innovations vector.
+% The authors have been contacted, fixed it, and report that the results
+% change only slightly at orderApp=3 to what they report in the paper. At
+% orderApp=2 all is correct and so the following part tests whether we get
+% the same model moments at the calibrated parameters (we do not optimize).
+% We compare it to the replication file RunGMM_Feedback_estim_RRA.m with the
+% following settings: orderApp=1|2, seOn=1, q_lag=10, weighting=1+1;
+% scaled=0; optimizer=0; estimator=1; momentSet=2;
+%
+% Output of the replication files for orderApp=1
+AndreasenEtAl.Q1 = 60275.3715;
+AndreasenEtAl.moments1 =[ % note that we reshuffeled to be compatible with our matched moments block
+ {[ 1]} {'Ex' } {'Gr_C '} {' ' } {'0.024388' } {'0.023726' }
+ {[ 2]} {'Ex' } {'Gr_I '} {' ' } {'0.031046' } {'0.027372' }
+ {[ 3]} {'Ex' } {'Infl ' } {' ' } {'0.03757' } {'0.041499' }
+ {[ 4]} {'Ex' } {'r1 ' } {' ' } {'0.056048' } {'0.077843' }
+ {[ 5]} {'Ex' } {'r40 ' } {' ' } {'0.069929' } {'0.077843' }
+ {[ 6]} {'Ex' } {'xhr40 '} {' ' } {'0.017237' } {'0' }
+ {[ 7]} {'Ex' } {'GoY '} {' ' } {'-1.5745' } {'-1.5746' }
+ {[ 8]} {'Ex' } {'hours '} {' ' } {'-0.043353' } {'-0.043299' }
+ {[ 9]} {'Exx' } {'Gr_C '} {'Gr_C '} {'0.0013159' } {'0.0012763' }
+ {[17]} {'Exx' } {'Gr_C '} {'Gr_I '} {'0.0021789' } {'0.0017759' }
+ {[18]} {'Exx' } {'Gr_C '} {'Infl ' } {'0.00067495' } {'0.00077354' }
+ {[19]} {'Exx' } {'Gr_C '} {'r1 ' } {'0.0011655' } {'0.0016538' }
+ {[20]} {'Exx' } {'Gr_C '} {'r40 ' } {'0.0015906' } {'0.0017949' }
+ {[21]} {'Exx' } {'Gr_C '} {'xhr40 '} {'0.0020911' } {'0.0017847' }
+ {[10]} {'Exx' } {'Gr_I '} {'Gr_I '} {'0.0089104' } {'0.0053424' }
+ {[22]} {'Exx' } {'Gr_I '} {'Infl ' } {'0.00063139' } {'0.00064897' }
+ {[23]} {'Exx' } {'Gr_I '} {'r1 ' } {'0.0011031' } {'0.0019533' }
+ {[24]} {'Exx' } {'Gr_I '} {'r40 ' } {'0.0018445' } {'0.0020602' }
+ {[25]} {'Exx' } {'Gr_I '} {'xhr40 '} {'0.00095556' } {'0.0064856' }
+ {[11]} {'Exx' } {'Infl ' } {'Infl ' } {'0.0020268' } {'0.0020922' }
+ {[26]} {'Exx' } {'Infl ' } {'r1 ' } {'0.0025263' } {'0.0036375' }
+ {[27]} {'Exx' } {'Infl ' } {'r40 ' } {'0.0029126' } {'0.0034139' }
+ {[28]} {'Exx' } {'Infl ' } {'xhr40 '} {'-0.00077101'} {'-0.0011665' }
+ {[12]} {'Exx' } {'r1 ' } {'r1 ' } {'0.0038708' } {'0.0066074' }
+ {[29]} {'Exx' } {'r1 ' } {'r40 ' } {'0.0044773' } {'0.0062959' }
+ {[30]} {'Exx' } {'r1 ' } {'xhr40 '} {'-0.00048202'} {'-0.00075499'}
+ {[13]} {'Exx' } {'r40 ' } {'r40 ' } {'0.0054664' } {'0.0061801' }
+ {[31]} {'Exx' } {'r40 ' } {'xhr40 '} {'0.00053864' } {'-0.00030456'}
+ {[14]} {'Exx' } {'xhr40 '} {'xhr40 '} {'0.053097' } {'0.012048' }
+ {[15]} {'Exx' } {'GoY '} {'GoY '} {'2.4863' } {'2.4872' }
+ {[16]} {'Exx' } {'hours '} {'hours '} {'0.0018799' } {'0.0018759' }
+ {[32]} {'Exx1'} {'Gr_C '} {'Gr_C '} {'0.00077917' } {'0.00080528' }
+ {[33]} {'Exx1'} {'Gr_I '} {'Gr_I '} {'0.0050104' } {'0.0017036' }
+ {[34]} {'Exx1'} {'Infl ' } {'Infl ' } {'0.0019503' } {'0.0020185' }
+ {[35]} {'Exx1'} {'r1 ' } {'r1 ' } {'0.0038509' } {'0.0065788' }
+ {[36]} {'Exx1'} {'r40 ' } {'r40 ' } {'0.0054699' } {'0.0061762' }
+ {[37]} {'Exx1'} {'xhr40 '} {'xhr40 '} {'-0.00098295'} {'-4.5519e-15'}
+ {[38]} {'Exx1'} {'GoY '} {'GoY '} {'2.4868' } {'2.4863' }
+ {[39]} {'Exx1'} {'hours '} {'hours '} {'0.0018799' } {'0.0018755' }
+];
+
+% Output of the replication files for orderApp=2
+AndreasenEtAl.Q2 = 140.8954;
+AndreasenEtAl.moments2 =[ % note that we reshuffeled to be compatible with our matched moments block
+ {[ 1]} {'Ex' } {'Gr_C '} {' ' } {'0.024388' } {'0.023726' }
+ {[ 2]} {'Ex' } {'Gr_I '} {' ' } {'0.031046' } {'0.027372' }
+ {[ 3]} {'Ex' } {'Infl ' } {' ' } {'0.03757' } {'0.034618' }
+ {[ 4]} {'Ex' } {'r1 ' } {' ' } {'0.056048' } {'0.056437' }
+ {[ 5]} {'Ex' } {'r40 ' } {' ' } {'0.069929' } {'0.07051' }
+ {[ 6]} {'Ex' } {'xhr40 '} {' ' } {'0.017237' } {'0.014242' }
+ {[ 7]} {'Ex' } {'GoY '} {' ' } {'-1.5745' } {'-1.574' }
+ {[ 8]} {'Ex' } {'hours '} {' ' } {'-0.043353' } {'-0.043351' }
+ {[ 9]} {'Exx' } {'Gr_C '} {'Gr_C '} {'0.0013159' } {'0.0012917' }
+ {[17]} {'Exx' } {'Gr_C '} {'Gr_I '} {'0.0021789' } {'0.0017862' }
+ {[18]} {'Exx' } {'Gr_C '} {'Infl ' } {'0.00067495' } {'0.00061078' }
+ {[19]} {'Exx' } {'Gr_C '} {'r1 ' } {'0.0011655' } {'0.0011494' }
+ {[20]} {'Exx' } {'Gr_C '} {'r40 ' } {'0.0015906' } {'0.0016149' }
+ {[21]} {'Exx' } {'Gr_C '} {'xhr40 '} {'0.0020911' } {'0.002203' }
+ {[10]} {'Exx' } {'Gr_I '} {'Gr_I '} {'0.0089104' } {'0.0054317' }
+ {[22]} {'Exx' } {'Gr_I '} {'Infl ' } {'0.00063139' } {'0.00045278' }
+ {[23]} {'Exx' } {'Gr_I '} {'r1 ' } {'0.0011031' } {'0.0013672' }
+ {[24]} {'Exx' } {'Gr_I '} {'r40 ' } {'0.0018445' } {'0.0018557' }
+ {[25]} {'Exx' } {'Gr_I '} {'xhr40 '} {'0.00095556' } {'0.0067742' }
+ {[11]} {'Exx' } {'Infl ' } {'Infl ' } {'0.0020268' } {'0.0016583' }
+ {[26]} {'Exx' } {'Infl ' } {'r1 ' } {'0.0025263' } {'0.0024521' }
+ {[27]} {'Exx' } {'Infl ' } {'r40 ' } {'0.0029126' } {'0.002705' }
+ {[28]} {'Exx' } {'Infl ' } {'xhr40 '} {'-0.00077101'} {'-0.00065007'}
+ {[12]} {'Exx' } {'r1 ' } {'r1 ' } {'0.0038708' } {'0.0038274' }
+ {[29]} {'Exx' } {'r1 ' } {'r40 ' } {'0.0044773' } {'0.004297' }
+ {[30]} {'Exx' } {'r1 ' } {'xhr40 '} {'-0.00048202'} {'6.3243e-05' }
+ {[13]} {'Exx' } {'r40 ' } {'r40 ' } {'0.0054664' } {'0.0051686' }
+ {[31]} {'Exx' } {'r40 ' } {'xhr40 '} {'0.00053864' } {'0.00066645' }
+ {[14]} {'Exx' } {'xhr40 '} {'xhr40 '} {'0.053097' } {'0.013543' }
+ {[15]} {'Exx' } {'GoY '} {'GoY '} {'2.4863' } {'2.4858' }
+ {[16]} {'Exx' } {'hours '} {'hours '} {'0.0018799' } {'0.0018804' }
+ {[32]} {'Exx1'} {'Gr_C '} {'Gr_C '} {'0.00077917' } {'0.00081772' }
+ {[33]} {'Exx1'} {'Gr_I '} {'Gr_I '} {'0.0050104' } {'0.0017106' }
+ {[34]} {'Exx1'} {'Infl ' } {'Infl ' } {'0.0019503' } {'0.0015835' }
+ {[35]} {'Exx1'} {'r1 ' } {'r1 ' } {'0.0038509' } {'0.0037985' }
+ {[36]} {'Exx1'} {'r40 ' } {'r40 ' } {'0.0054699' } {'0.0051642' }
+ {[37]} {'Exx1'} {'xhr40 '} {'xhr40 '} {'-0.00098295'} {'0.00020285' }
+ {[38]} {'Exx1'} {'GoY '} {'GoY '} {'2.4868' } {'2.4848' }
+ {[39]} {'Exx1'} {'hours '} {'hours '} {'0.0018799' } {'0.0018799' }
+];
+
+@#for orderApp in 1:2
+
+method_of_moments(
+ mom_method = GMM % method of moments method; possible values: GMM|SMM
+ , datafile = 'AFVRR_data.mat' % name of filename with data
+ , bartlett_kernel_lag = 10 % bandwith in optimal weighting matrix
+ , order = @{orderApp} % order of Taylor approximation in perturbation
+ , pruning % use pruned state space system at higher-order
+ % , verbose % display and store intermediate estimation results
+ , weighting_matrix = ['DIAGONAL'] % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
+ % , TeX % print TeX tables and graphics
+ % Optimization options that can be set by the user in the mod file, otherwise default values are provided
+ %, huge_number=1D10 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+ , mode_compute = 0 % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
+ , optim = ('TolFun', 1e-6
+ ,'TolX', 1e-6
+ ,'MaxIter', 3000
+ ,'MaxFunEvals', 1D6
+ ,'UseParallel' , 1
+ %,'Jacobian' , '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_standard_errors
+ , se_tolx=1e-10
+);
+
+% Check results
+
+fprintf('****************************************************************\n')
+fprintf('Compare Results for perturbation order @{orderApp}\n')
+fprintf('****************************************************************\n')
+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;
+
+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))])
+
+if norm(dev_modelmoments)> 1e-4
+ warning('Something wrong in the computation of moments at order @{orderApp}')
+end
+
+@#endfor
+
+%--------------------------------------------------------------------------
+% Replicate estimation at orderApp=3
+%--------------------------------------------------------------------------
+@#ifdef DoEstimation
+method_of_moments(
+ mom_method = GMM % method of moments method; possible values: GMM|SMM
+ , datafile = 'AFVRR_data.mat' % name of filename with data
+ , bartlett_kernel_lag = 10 % bandwith in optimal weighting matrix
+ , order = 3 % order of Taylor approximation in perturbation
+ , pruning % use pruned state space system at higher-order
+ % , verbose % display and store intermediate estimation results
+ , weighting_matrix = ['DIAGONAL', 'Optimal'] % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
+ % , TeX % print TeX tables and graphics
+ % Optimization options that can be set by the user in the mod file, otherwise default values are provided
+ %, huge_number=1D10 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+ , mode_compute = 13 % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
+ , additional_optimizer_steps = [13]
+ , optim = ('TolFun', 1e-6
+ ,'TolX', 1e-6
+ ,'MaxIter', 3000
+ ,'MaxFunEvals', 1D6
+ ,'UseParallel' , 1
+ %,'Jacobian' , '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_standard_errors
+ , se_tolx=1e-10
+);
+@#endif
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/AFVRR/AFVRR_common.inc b/tests/estimation/method_of_moments/AFVRR/AFVRR_common.inc
new file mode 100644
index 000000000..76aea9e0b
--- /dev/null
+++ b/tests/estimation/method_of_moments/AFVRR/AFVRR_common.inc
@@ -0,0 +1,540 @@
+% DSGE model based on replication files of
+% Andreasen, Fernandez-Villaverde, Rubio-Ramirez (2018), The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications, Review of Economic Studies, 85, p. 1-49
+% Original code by Martin M. Andreasen, Jan 2016
+% Adapted for Dynare by Willi Mutschler (@wmutschl, willi@mutschler.eu), Jan 2021
+% =========================================================================
+% 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 .
+% =========================================================================
+
+%--------------------------------------------------------------------------
+% Variable declaration
+%--------------------------------------------------------------------------
+var
+ln_k
+ln_s
+ln_a
+ln_g
+ln_d
+
+ln_c
+ln_r
+ln_pai
+ln_h
+ln_q
+ln_evf
+ln_iv
+ln_x2
+ln_la
+ln_goy
+ln_Esdf
+
+xhr20
+xhr40
+Exhr
+
+@#for i in 1:40
+ln_p@{i}
+@#endfor
+
+Obs_Gr_C
+Obs_Gr_I
+Obs_Infl
+Obs_r1
+Obs_r40
+Obs_xhr40
+Obs_GoY
+Obs_hours
+;
+
+predetermined_variables ln_k ln_s;
+
+varobs Obs_Gr_C Obs_Gr_I Obs_Infl Obs_r1 Obs_r40 Obs_xhr40 Obs_GoY Obs_hours;
+
+%--------------------------------------------------------------------------
+% Exogenous shocks
+%--------------------------------------------------------------------------
+varexo
+eps_a
+eps_d
+eps_g
+;
+
+%--------------------------------------------------------------------------
+% Parameter declaration
+%--------------------------------------------------------------------------
+parameters
+BETTA
+B
+INHABIT
+H
+PHI1
+PHI2
+RRA
+PHI4
+KAPAone
+KAPAtwo
+DELTA
+THETA
+ETA
+ALFA
+CHI
+RHOR
+BETTAPAI
+BETTAY
+MYYPS
+MYZ
+RHOA
+%STDA
+RHOG
+%STDG
+RHOD
+%STDD
+CONSxhr40
+BETTAxhr
+BETTAxhr40
+CONSxhr20
+PAI
+GAMA
+GoY
+
+%auxiliary
+PHIzero
+AA
+PHI3
+negVf
+;
+
+
+%--------------------------------------------------------------------------
+% Model equations
+%--------------------------------------------------------------------------
+% Based on DSGE_model_NegVf_yieldCurve.m and DSGE_model_PosVf_yieldCurve.m
+% Note that we include an auxiliary parameter negVf to distinguish whether
+% the steady state value function is positive (negVf=0) or negative (negVf=1).
+% This parameter is endogenously determined in the steady_state_model block.
+
+model;
+%--------------------------------------------------------------------------
+% Auxiliary expressions
+%--------------------------------------------------------------------------
+% do exp transform such that variables are logged variables
+@#for var in [ "k", "s", "c", "r", "a", "g", "d", "pai", "h", "q", "evf", "iv", "x2", "la", "goy", "Esdf" ]
+#@{var}_ba1 = exp(ln_@{var}(-1));
+#@{var}_cu = exp(ln_@{var});
+#@{var}_cup = exp(ln_@{var}(+1));
+@#endfor
+@#for i in 1:40
+#p@{i}_cu = exp(ln_p@{i});
+#p@{i}_cup = exp(ln_p@{i}(+1));
+@#endfor
+% these variables are not transformed
+#xhr20_cu = xhr20;
+#xhr20_cup = xhr20(+1);
+#xhr40_cu = xhr40;
+#xhr40_cup = xhr40(+1);
+#Exhr_cu = Exhr;
+#Exhr_cup = Exhr(+1);
+
+% auxiliary steady state variables
+#K = exp(steady_state(ln_k));
+#IV = exp(steady_state(ln_iv));
+#C = exp(steady_state(ln_c));
+#Y = (C + IV)/(1-GoY);
+#R = exp(steady_state(ln_r));
+#G = Y-C-IV;
+
+#removeMeanXhr = 1;
+
+% The atemporal relations if possible
+% No stochastic trend in investment specific shocks
+#myyps_cu = MYYPS;
+#myyps_cup = MYYPS;
+
+% No stochastic trend in non-stationary technology shocks
+#myz_cu = MYZ;
+#myz_cup = MYZ;
+
+% Defining myzstar
+#MYZSTAR = MYYPS^(THETA/(1-THETA))*MYZ;
+#myzstar_cu = myyps_cu ^(THETA/(1-THETA))*myz_cu;
+#myzstar_cup= myyps_cup^(THETA/(1-THETA))*myz_cup;
+
+% The expression for the value function (only valid for deterministic trends!)
+% Note that we make use of auxiliary parameter negVf to switch signs
+#mvf_cup = -negVf*(d_cup/(1-PHI2)*((c_cup-B*c_cu*MYZSTAR^-1)^(1-PHI2)-1) + d_cup*PHIzero/(1-PHI1)*(1-h_cup)^(1-PHI1) - negVf* BETTA*MYZSTAR^((1-PHI4)*(1-PHI2))*AA*evf_cup^(1/(1-PHI3)));
+
+% The growth rate in lambda
+#myla_cup = (la_cup/la_cu)*(AA*evf_cu^(1/(1-PHI3))/mvf_cup)^PHI3*myzstar_cup^(-PHI2*(1-PHI4)-PHI4);
+
+% The relation between the optimal price for the firms and the pris and inflation
+%ptil_cu = ((1-ALFA*(pai_ba1^CHI/pai_cu )^(1-ETA))/(1-ALFA))^(1/(1-ETA));
+%ptil_cup = ((1-ALFA*(pai_cu ^CHI/pai_cup)^(1-ETA))/(1-ALFA))^(1/(1-ETA));
+#ptil_cu = ((1-ALFA*(1/pai_cu )^(1-ETA))/(1-ALFA))^(1/(1-ETA));
+#ptil_cup = ((1-ALFA*(1/pai_cup)^(1-ETA))/(1-ALFA))^(1/(1-ETA));
+
+% From the households' FOC for labor
+#w_cu = d_cu*PHIzero*(1-h_cu )^(-PHI1)/la_cu;
+#w_cup = d_cu*PHIzero*(1-h_cup)^(-PHI1)/la_cup;
+% Shouldn't w_cup include d_cup? Let's stick to the original (wrong) code in the replication files as results don't change dramatically... [@wmutschl]
+
+% The firms' FOC for labor
+#mc_cu = w_cu /((1-THETA)*a_cu *myyps_cu ^(-THETA/(1-THETA))*myz_cu ^-THETA *k_cu ^THETA*h_cu ^(-THETA));
+#mc_cup = w_cup/((1-THETA)*a_cup*myyps_cup^(-THETA/(1-THETA))*myz_cup^-THETA *k_cup^THETA*h_cup^(-THETA));
+
+% The firms' FOC for capital
+#rk_cu = mc_cu *THETA* a_cu *myyps_cu *myz_cu ^(1-THETA)*k_cu ^(THETA-1)*h_cu ^(1-THETA);
+#rk_cup = mc_cup*THETA* a_cup*myyps_cup*myz_cup^(1-THETA)*k_cup^(THETA-1)*h_cup^(1-THETA);
+
+% The income identity
+#y_cu = c_cu + iv_cu + g_cu;
+
+%--------------------------------------------------------------------------
+% Actual model equations
+%--------------------------------------------------------------------------
+
+[name='Expected value of the value function']
+0 = -evf_cu + (mvf_cup/AA)^(1-PHI3);
+
+[name='Households FOC for capital']
+0 = -q_cu+BETTA*myla_cup/myyps_cup*(rk_cup+q_cup*(1-DELTA) -q_cup*KAPAtwo/2*(iv_cup/k_cup*myyps_cup*myzstar_cup - IV/K*MYYPS*MYZSTAR)^2 +q_cup*KAPAtwo*(iv_cup/k_cup*myyps_cup*myzstar_cup - IV/K*MYYPS*MYZSTAR)*iv_cup/k_cup*myyps_cup*myzstar_cup);
+
+[name='Households FOC for investments']
+0 = -1+q_cu*(1-KAPAone/2*(iv_cu/IV-1)^2-iv_cu/IV*KAPAone*(iv_cu/IV-1)-KAPAtwo*(iv_cu/k_cu*myyps_cu*myzstar_cu - IV/K*MYYPS*MYZSTAR));
+
+[name='Euler equation for consumption']
+0 = -1+BETTA*r_cu*exp(CONSxhr40*xhr40_cu + CONSxhr20*xhr20_cu)*myla_cup/pai_cup;
+
+[name='Households FOC for consumption']
+0 = -la_cu + d_cu*(c_cu -B*c_ba1*myzstar_cu^-1)^(-PHI2) -INHABIT*B*BETTA*d_cup*(AA*evf_cu^(1/(1-PHI3))/mvf_cup)^PHI3*(c_cup -B*c_cu*myzstar_cup^-1)^(-PHI2)*myzstar_cup^(-PHI2*(1-PHI4)-PHI4);
+
+[name='Nonlinear pricing, relation for x1 = (ETA-1)/ETA*x2']
+0= -(ETA-1)/ETA*x2_cu+y_cu*mc_cu*ptil_cu^(-ETA-1) +ALFA*BETTA*myla_cup*(ptil_cu/ptil_cup)^(-ETA-1)*(1/pai_cup)^(-ETA)*(ETA-1)/ETA*x2_cup*myzstar_cup;
+
+[name='Nonlinear pricing, relation for x2']
+0=-x2_cu+y_cu*ptil_cu^-ETA +ALFA*BETTA*myla_cup*(ptil_cu/ptil_cup)^(-ETA)*(1/pai_cup)^(1-ETA)*x2_cup*myzstar_cup;
+
+[name='Nonlinear pricing, relation for s']
+0= -s_cup+(1-ALFA)*ptil_cu^(-ETA)+ALFA*(pai_cu/1)^ETA*s_cu;
+
+[name='Interest rate rule']
+0 = -log(r_cu/R)+RHOR*log(r_ba1/R)+(1-RHOR)*(BETTAPAI*log(pai_cu/PAI)+BETTAY*log(y_cu/Y) + BETTAxhr*(BETTAxhr40*xhr40_cu - removeMeanXhr*Exhr_cu));
+
+[name='Production function']
+0 = -y_cu*s_cup + a_cu *(k_cu *myyps_cu ^(-1/(1-THETA))*myz_cu ^-1)^THETA*h_cu ^(1-THETA);
+
+[name='Relation for physical capital stock']
+0= -k_cup + (1-DELTA)*k_cu*(myyps_cu*myzstar_cu)^-1 + iv_cu - iv_cu*KAPAone/2*(iv_cu/IV-1)^2 - k_cu*(myyps_cu*myzstar_cu)^-1*KAPAtwo/2*(iv_cu/k_cu*myyps_cu*myzstar_cu - IV/K*MYYPS*MYZSTAR)^2;
+
+[name='Goverment spending over output']
+0=-goy_cu + g_cu/y_cu;
+
+[name='The yield curve: p1']
+0= -p1_cu + 1/r_cu;
+
+@#for i in 2:40
+[name='The yield curve: p@{i}']
+0= -p@{i}_cu + BETTA*myla_cup/pai_cup*p@{i-1}_cup;
+@#endfor
+
+[name='Stochastic discount factor']
+0= -Esdf_cu+ BETTA*myla_cup/pai_cup;
+
+[name='Expected 5 year excess holding period return']
+0= -xhr20_cu+ log(p19_cup) - log(p20_cu) - log(r_cu);
+
+[name='Expected 10 year excess holding period return']
+0= -xhr40_cu+ log(p39_cup) - log(p40_cu) - log(r_cu);
+
+[name='Mean of expected excess holding period return in Taylor rule']
+0= -Exhr_cu + (1-GAMA)*(BETTAxhr40*xhr40_cu) + GAMA*Exhr_cup;
+
+[name='Exogenous process for productivity']
+0 = -log(a_cu)+RHOA*log(a_ba1) + eps_a;
+
+[name='Exogenous process for government spending']
+0 = -log(g_cu/G)+RHOG*log(g_ba1/G) + eps_g;
+
+[name='Exogenous process for discount factor shifter']
+0 = -log(d_cu)+RHOD*log(d_ba1) + eps_d;
+
+[name='Observable annualized consumption growth']
+Obs_Gr_C = 4*( ln_c -ln_c(-1) + log(MYZSTAR));
+
+[name='Observable annualized investment growth']
+Obs_Gr_I = 4*( ln_iv - ln_iv(-1) + log(MYZSTAR)+log(MYYPS));
+
+[name='Observable annualized inflation']
+Obs_Infl = 4*( ln_pai);
+
+[name='Observable annualized one-quarter nominal yield']
+Obs_r1 = 4*( ln_r);
+
+[name='Observable annualized 10-year nominal yield']
+Obs_r40 = 4*( -1/40*ln_p40);
+
+[name='Observable annualized 10-year ex post excess holding period return']
+Obs_xhr40 = 4*( ln_p39 - ln_p40(-1) - ln_r(-1) );
+
+[name='Observable annualized log ratio of government spending to GDP']
+Obs_GoY = 4*( 1/4*ln_goy);
+
+[name='Observable annualized log of hours']
+Obs_hours = 4*( 1/100*ln_h);
+end;
+
+
+%--------------------------------------------------------------------------
+% Steady State Computations
+%--------------------------------------------------------------------------
+% Based on DSGE_model_yieldCurve_ss.m, getPHI3.m, ObjectGMM.m
+% Note that we include an auxiliary parameter negVf to distinguish whether
+% the steady state value function is positive (negVf=0) or negative (negVf=1).
+% This parameter is endogenously determined in the steady_state_model block.
+
+
+steady_state_model;
+
+% The growth rate in the firms' fixed costs
+MYZSTARBAR = MYYPS^(THETA/(1-THETA))*MYZ;
+
+% The growth rate for lampda
+MYLABAR = MYZSTARBAR^(-PHI2*(1-PHI4)-PHI4);
+
+% The relative optimal price for firms
+PTILBAR = ((1-ALFA*PAI^((CHI-1)*(1-ETA)))/(1-ALFA))^(1/(1-ETA));
+
+% The state variable s for distortions between output and produktion
+SBAR = ((1-ALFA)*PTILBAR^(-ETA))/(1-ALFA*PAI^((1-CHI)*ETA));
+
+% The 1-period interest rate
+RBAR = PAI/(BETTA*MYLABAR);
+
+% The market price of capital
+QBAR = 1;
+
+% The real price of renting capital
+RKBAR = QBAR*(MYYPS/(BETTA*MYLABAR)-(1-DELTA));
+
+% The marginal costs in the firms
+MCBAR = (1-ALFA*BETTA*MYLABAR*PAI^((1-CHI)*ETA)*MYZSTARBAR)*(ETA-1)/ETA*PTILBAR/(1-ALFA*BETTA*MYLABAR*PAI^((CHI-1)*(1-ETA))*MYZSTARBAR);
+
+% The capital stock
+KBAR = H*(RKBAR/(MCBAR*THETA*MYYPS*MYZ^(1-THETA)))^(1/(THETA-1));
+
+% The wage level
+WBAR = MCBAR*(1-THETA)*MYYPS^(-THETA/(1-THETA))*MYZ^-THETA*(KBAR/H)^THETA;
+
+% The level of investment
+IVBAR = KBAR - (1-DELTA)*KBAR*MYYPS^(-1/(1-THETA))*MYZ^-1;
+
+% The consumption level
+CBAR = ((1-GoY)*(KBAR*MYYPS^(-1/(1-THETA))*MYZ^-1)^THETA*H^(1-THETA))/SBAR-IVBAR;
+
+% The output level
+YBAR = (CBAR + IVBAR)/(1-GoY);
+
+% The value of lambda
+LABAR = (CBAR-B*CBAR*MYZSTARBAR^-1)^-PHI2 - INHABIT*B*BETTA*(CBAR-B*CBAR*MYZSTARBAR^-1)^-PHI2*MYZSTARBAR^(-PHI2*(1-PHI4)-PHI4);
+
+% The value of PHIzero
+PHIzero = LABAR*WBAR*(1-H)^PHI1;
+
+% The level of the value function
+VFBAR = 1/(1-BETTA*MYZSTARBAR^((1-PHI4)*(1-PHI2)))*(1/(1-PHI2)*((CBAR-B*CBAR*MYZSTARBAR^-1)^(1-PHI2)-1)+PHIzero/(1-PHI1)*(1-H)^(1-PHI1));
+UBAR = 1/(1-PHI2)*((CBAR-B*CBAR*MYZSTARBAR^-1)^(1-PHI2)-1)+PHIzero/(1-PHI1)*(1-H)^(1-PHI1);
+[AA, EVFBAR, PHI3, negVf, info]= AFVRR_steady_helper(VFBAR,RBAR,IVBAR,CBAR,KBAR,LABAR,QBAR,YBAR, BETTA,B,PAI,H,PHIzero,PHI1,PHI2,THETA,MYYPS,MYZ,INHABIT,RRA,CONSxhr40);
+% The value of X2
+X2BAR = YBAR*PTILBAR^(-ETA)/(1-BETTA*ALFA*MYLABAR*PAI^((CHI-1)*(1-ETA))*MYZSTARBAR);
+
+% Government spending
+GBAR = GoY*YBAR;
+%**************************************************************************
+
+% map into model variables
+ln_k = log(KBAR);
+ln_s = log(SBAR);
+ln_c_ba1 = log(CBAR);
+ln_r_ba1 = log(RBAR);
+ln_a = log(1);
+ln_g = log(GBAR);
+ln_d = log(1);
+
+ln_c = log(CBAR);
+ln_r = log(RBAR);
+ln_pai = log(PAI);
+ln_h = log(H);
+ln_q = log(QBAR);
+ln_evf = log(EVFBAR);
+ln_iv = log(IVBAR);
+ln_x2 = log(X2BAR);
+ln_la = log(LABAR);
+ln_goy = log(GoY);
+ln_Esdf = log(1/RBAR);
+xhr20 = 0;
+xhr40 = 0;
+Exhr = 0;
+% The yield curve
+ln_p1 = log((1/RBAR)^1);
+ln_p2 = log((1/RBAR)^2);
+ln_p3 = log((1/RBAR)^3);
+ln_p4 = log((1/RBAR)^4);
+ln_p5 = log((1/RBAR)^5);
+ln_p6 = log((1/RBAR)^6);
+ln_p7 = log((1/RBAR)^7);
+ln_p8 = log((1/RBAR)^8);
+ln_p9 = log((1/RBAR)^9);
+ln_p10 = log((1/RBAR)^10);
+ln_p11 = log((1/RBAR)^11);
+ln_p12 = log((1/RBAR)^12);
+ln_p13 = log((1/RBAR)^13);
+ln_p14 = log((1/RBAR)^14);
+ln_p15 = log((1/RBAR)^15);
+ln_p16 = log((1/RBAR)^16);
+ln_p17 = log((1/RBAR)^17);
+ln_p18 = log((1/RBAR)^18);
+ln_p19 = log((1/RBAR)^19);
+ln_p20 = log((1/RBAR)^20);
+ln_p21 = log((1/RBAR)^21);
+ln_p22 = log((1/RBAR)^22);
+ln_p23 = log((1/RBAR)^23);
+ln_p24 = log((1/RBAR)^24);
+ln_p25 = log((1/RBAR)^25);
+ln_p26 = log((1/RBAR)^26);
+ln_p27 = log((1/RBAR)^27);
+ln_p28 = log((1/RBAR)^28);
+ln_p29 = log((1/RBAR)^29);
+ln_p30 = log((1/RBAR)^30);
+ln_p31 = log((1/RBAR)^31);
+ln_p32 = log((1/RBAR)^32);
+ln_p33 = log((1/RBAR)^33);
+ln_p34 = log((1/RBAR)^34);
+ln_p35 = log((1/RBAR)^35);
+ln_p36 = log((1/RBAR)^36);
+ln_p37 = log((1/RBAR)^37);
+ln_p38 = log((1/RBAR)^38);
+ln_p39 = log((1/RBAR)^39);
+ln_p40 = log((1/RBAR)^40);
+
+Obs_Gr_C = 4*( log(MYZSTARBAR) );
+Obs_Gr_I = 4*( log(MYZSTARBAR)+log(MYYPS) );
+Obs_Infl = 4*( ln_pai );
+Obs_r1 = 4*( ln_r );
+Obs_r40 = 4*( -1/40*ln_p40 );
+Obs_xhr40 = 4*( xhr40 );
+Obs_GoY = 4*( 1/4*ln_goy );
+Obs_hours = 4*( 1/100*ln_h );
+end;
+
+%--------------------------------------------------------------------------
+% Declare moments to use in estimation
+%--------------------------------------------------------------------------
+% These are the moments used in the paper; corresponds to momentSet=2 in the replication files
+
+matched_moments;
+%mean
+Obs_Gr_C;
+Obs_Gr_I;
+Obs_Infl;
+Obs_r1;
+Obs_r40;
+Obs_xhr40;
+Obs_GoY;
+Obs_hours;
+
+% all variances
+Obs_Gr_C*Obs_Gr_C;
+Obs_Gr_I*Obs_Gr_I;
+Obs_Infl*Obs_Infl;
+Obs_r1*Obs_r1;
+Obs_r40*Obs_r40;
+Obs_xhr40*Obs_xhr40;
+Obs_GoY*Obs_GoY;
+Obs_hours*Obs_hours;
+
+% covariance excluding GoY and hours
+Obs_Gr_C*Obs_Gr_I;
+Obs_Gr_C*Obs_Infl;
+Obs_Gr_C*Obs_r1;
+Obs_Gr_C*Obs_r40;
+Obs_Gr_C*Obs_xhr40;
+%Obs_Gr_C*Obs_GoY;
+%Obs_Gr_C*Obs_hours;
+
+Obs_Gr_I*Obs_Infl;
+Obs_Gr_I*Obs_r1;
+Obs_Gr_I*Obs_r40;
+Obs_Gr_I*Obs_xhr40;
+%Obs_Gr_I*Obs_GoY;
+%Obs_Gr_I*Obs_hours;
+
+Obs_Infl*Obs_r1;
+Obs_Infl*Obs_r40;
+Obs_Infl*Obs_xhr40;
+%Obs_Infl*Obs_GoY;
+%Obs_Infl*Obs_hours;
+
+Obs_r1*Obs_r40;
+Obs_r1*Obs_xhr40;
+%Obs_r1*Obs_GoY;
+%Obs_r1*Obs_hours;
+
+Obs_r40*Obs_xhr40;
+%Obs_r40*Obs_GoY;
+%Obs_r40*Obs_hours;
+
+%Obs_xhr40*Obs_GoY;
+%Obs_xhr40*Obs_hours;
+
+%Obs_GoY*Obs_hours;
+
+%first autocovariance
+Obs_Gr_C*Obs_Gr_C(-1);
+Obs_Gr_I*Obs_Gr_I(-1);
+Obs_Infl*Obs_Infl(-1);
+Obs_r1*Obs_r1(-1);
+Obs_r40*Obs_r40(-1);
+Obs_xhr40*Obs_xhr40(-1);
+Obs_GoY*Obs_GoY(-1);
+Obs_hours*Obs_hours(-1);
+end;
+
+%--------------------------------------------------------------------------
+% Create Data
+%--------------------------------------------------------------------------
+@#ifdef CreateData
+verbatim;
+% From 1961Q3 to 2007Q4
+DataUS = xlsread('Data_PruningPaper_v5.xlsx','Data_used','E3:M188');
+% ANNUALIZED (except for hours and GoY)
+% 1 2 3 4 5 6 7 8 9
+% Lables: Date Gr_C Gr_I GoY hours Infl_C r1 r40 xhr40
+%label_data = {'Gr_C ', 'Gr_I ','Infl ', 'r1 ', 'r40 ', 'xhr40 ','GoY ', 'hours '};
+%DataUS = [DataUS(:,2:3) DataUS(:,6:8) DataUS(:,9) log(DataUS(:,4)) 4*log(DataUS(:,5))/100];
+Obs_Gr_C = DataUS(:,2);
+Obs_Gr_I = DataUS(:,3);
+Obs_Infl = DataUS(:,6);
+Obs_r1 = DataUS(:,7);
+Obs_r40 = DataUS(:,8);
+Obs_xhr40 = DataUS(:,9);
+Obs_GoY = log(DataUS(:,4));
+Obs_hours = 4*log(DataUS(:,5))/100;
+
+save('AFVRR_data.mat','Obs_Gr_C','Obs_Gr_I','Obs_Infl','Obs_r1','Obs_r40','Obs_xhr40','Obs_GoY','Obs_hours');
+pause(1);
+end;
+@#endif
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/AFVRR/AFVRR_data.mat b/tests/estimation/method_of_moments/AFVRR/AFVRR_data.mat
new file mode 100644
index 000000000..f606b2109
Binary files /dev/null and b/tests/estimation/method_of_moments/AFVRR/AFVRR_data.mat differ
diff --git a/tests/estimation/method_of_moments/AFVRR/AFVRR_steady_helper.m b/tests/estimation/method_of_moments/AFVRR/AFVRR_steady_helper.m
new file mode 100644
index 000000000..b8289d484
--- /dev/null
+++ b/tests/estimation/method_of_moments/AFVRR/AFVRR_steady_helper.m
@@ -0,0 +1,80 @@
+% DSGE model based on replication files of
+% Andreasen, Fernandez-Villaverde, Rubio-Ramirez (2018), The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications, Review of Economic Studies, 85, p. 1-49
+% Adapted for Dynare by Willi Mutschler (@wmutschl, willi@mutschler.eu), Jan 2021
+% =========================================================================
+% 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 .
+% =========================================================================
+
+% This is a helper function to compute steady state values and endogenous parameters
+% Based on DSGE_model_yieldCurve_ss.m, getPHI3.m, ObjectGMM.m
+
+function [AA, EVFBAR, PHI3, negVf, info]= AFVRR_steady_helper(VFBAR,RBAR,IVBAR,CBAR,KBAR,LABAR,QBAR,YBAR, BETTA,B,PAI,H,PHIzero,PHI1,PHI2,THETA,MYYPS,MYZ,INHABIT,RRA,CONSxhr40)
+% We get nice values of EVF by setting AA app. equal to VF.
+% The value of the expected value function raised to the power 1-PHI3
+% Also we check bounds on other variables
+% % Adding PHI3 to params. Note that PHI3 only affects the value function in
+% % steady state, hence the value we assign to PHI3 is irrelevant
+% PHI3 = -100;
+
+info=0;
+AA = NaN;
+EVFBAR = NaN;
+PHI3 = NaN;
+negVf = NaN;
+
+MYZSTAR = MYYPS^(THETA/(1-THETA))*MYZ;
+% The wage level
+WBAR = PHIzero*(1-H)^(-PHI1)/LABAR;
+RRAc = RRA;
+if INHABIT == 1
+ PHI3 = (RRAc - PHI2/((1-B*MYZSTAR^-1)/(1-BETTA*B)+PHI2/PHI1*WBAR*(1-H)/CBAR))/((1-PHI2)/((1-B*MYZSTAR^-1)/(1-BETTA*B)-(CBAR-B*CBAR*MYZSTAR^-1)^PHI2/((1-BETTA*B)*CBAR)+WBAR*(1-H)/CBAR*(1-PHI2)/(1-PHI1)));
+else
+ PHI3 = (RRAc - PHI2/(1-B*MYZSTAR^-1+PHI2/PHI1*WBAR*(1-H)/CBAR))/((1-PHI2)/(1-B*MYZSTAR^-1-(CBAR-B*CBAR*MYZSTAR^-1)^PHI2/((1-BETTA*B)*CBAR)+WBAR*(1-H)/CBAR*(1-PHI2)/(1-PHI1)));
+end
+if abs(PHI3) > 30000
+ disp('abs of PHI3 exceeds 30000')
+ info=1;
+ return
+end
+
+if CONSxhr40 > 1
+ info=1;
+ return
+end
+
+
+if VFBAR < 0
+ AA = -VFBAR;
+ EVFBAR = (-VFBAR/AA)^(1-PHI3);
+ negVf = 1;
+else
+ AA = VFBAR;
+ EVFBAR = (VFBAR/AA)^(1-PHI3);
+ negVf = -1;
+ disp('Positive Value Function');
+end
+
+
+if RBAR < 1 || IVBAR < 0 || CBAR < 0 || KBAR < 0 || PAI < 1 || H < 0 || H > 1 || QBAR < 0 || YBAR < 0
+ info = 1;
+end
+
+end
+
+
+
diff --git a/tests/estimation/method_of_moments/AnScho_MoM.mod b/tests/estimation/method_of_moments/AnScho/AnScho_MoM.mod
similarity index 73%
rename from tests/estimation/method_of_moments/AnScho_MoM.mod
rename to tests/estimation/method_of_moments/AnScho/AnScho_MoM.mod
index 83cec9136..3d0e718c9 100644
--- a/tests/estimation/method_of_moments/AnScho_MoM.mod
+++ b/tests/estimation/method_of_moments/AnScho/AnScho_MoM.mod
@@ -1,8 +1,8 @@
-% DSGE model used in replication files of
+% DSGE model used in replication files of
% An, Sungbae and Schorfheide, Frank, (2007), Bayesian Analysis of DSGE Models, Econometric Reviews, 26, issue 2-4, p. 113-172.
% Adapted by Willi Mutschler (@wmutschl, willi@mutschler.eu)
% =========================================================================
-% Copyright (C) 2020 Dynare Team
+% Copyright (C) 2020-2021 Dynare Team
%
% This file is part of Dynare.
%
@@ -203,28 +203,33 @@ end
@#for mommethod in ["GMM", "SMM"]
method_of_moments(
- % Necessery options
- mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM
- , datafile = 'AnScho_MoM_data_@{orderApp}.mat' % name of filename with data
+ % Necessery options
+ mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM
+ , datafile = 'AnScho_MoM_data_@{orderApp}.mat' % name of filename with data
- % Options for both GMM and SMM
+ % Options for both GMM and SMM
% , bartlett_kernel_lag = 20 % bandwith in optimal weighting matrix
, order = @{orderApp} % order of Taylor approximation in perturbation
- % , penalized_estimator % use penalized optimization
+ % , penalized_estimator % include deviation from prior mean as additional moment restriction and use prior precision as weight
, pruning % use pruned state space system at higher-order
% , verbose % display and store intermediate estimation results
- , weighting_matrix = ['optimal'] % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
- , additional_optimizer_steps = [4] % vector of numbers for the iterations in the 2-step feasible method of moments
- % , prefilter=0 % demean each data series by its empirical mean and use centered moments
- %
- % Options for SMM
+ , 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']
+ %, weighting_matrix_scaling_factor=1 % scaling of weighting matrix in objective function
+ , se_tolx=1e-6 % step size for numerical computation of standard errors
+
+ % Options for SMM
+ % , burnin=500 % number of periods dropped at beginning of simulation
% , bounded_shock_support % trim shocks in simulation to +- 2 stdev
- % , drop = 500 % number of periods dropped at beginning of simulation
% , seed = 24051986 % seed used in simulations
% , simulation_multiple = 5 % multiple of the data length used for simulation
- %
- % General options
- %, dirname = 'MM' % directory in which to store estimation output
+
+ % Options for GMM
+ @#if mommethod == "GMM"
+ , analytic_standard_errors % compute standard errors using analytical derivatives
+ @#endif
+
+ % General options
+ % , dirname = 'MM' % directory in which to store estimation output
% , graph_format = EPS % specify the file format(s) for graphs saved to disk
% , nodisplay % do not display the graphs, but still save them to disk
% , nograph % do not create graphs (which implies that they are not saved to the disk nor displayed)
@@ -232,44 +237,50 @@ end
% , plot_priors = 1 % control plotting of priors
% , prior_trunc = 1e-10 % probability of extreme values of the prior density that is ignored when computing bounds for the parameters
% , TeX % print TeX tables and graphics
- %
- % Data and model options
- %, first_obs = 501 % number of first observation
- % , logdata % if loglinear is set, this option is necessary if the user provides data already in logs, otherwise the log transformation will be applied twice (this may result in complex data)
- % , loglinear % computes a log-linear approximation of the model instead of a linear approximation
- , nobs = 250 % number of observations
- % , xls_sheet = willi % name of sheet with data in Excel
+
+ % Data and model options
+ % , first_obs = 501 % number of first observation
+ % , logdata % if data is already in logs
+ , nobs = 250 % number of observations
+ % , prefilter=0 % demean each data series by its empirical mean and use centered moments
+
+ % , xls_sheet = data % name/number of sheet with data in Excel
% , xls_range = B2:D200 % range of data in Excel sheet
- %
- % Optimization options that can be set by the user in the mod file, otherwise default values are provided
- % , analytic_derivation % uses analytic derivatives to compute standard errors for GMM
- %, huge_number=1D10 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
- , mode_compute = @{optimizer} % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
- %, optim = ('TolFun', 1e-5
- % ,'TolX', 1e-6
- % ) % 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
- % , tolf = 1e-5 % convergence criterion on function value for numerical differentiation
- % , tolx = 1e-6 % convergence criterion on funciton input for numerical differentiation
- %
- % % Numerical algorithms options
+
+ % Optimization options that can be set by the user in the mod file, otherwise default values are provided
+ % , huge_number=1e7 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+ , mode_compute = @{optimizer} % specifies the optimizer for minimization of moments distance
+ , additional_optimizer_steps = [1] % vector of additional mode-finders run after mode_compute
+ % optim: a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute, some exemplary common options:
+ , optim = ('TolFun' , 1e-6 % termination tolerance on the function value, a positive scalar
+ ,'TolX' , 1e-6 % termination tolerance on x, a positive scalar
+ ,'MaxIter' , 3000 % maximum number of iterations allowed, a positive integer
+ ,'MaxFunEvals' , 1D6 % maximum number of function evaluations allowed, a positive integer
+ % ,'UseParallel' , 1 % when true (and supported by optimizer) solver estimates gradients in parallel (using Matlab/Octave's parallel toolbox)
+ % ,'Jacobian' , 'off' % when 'off' gradient-based solvers approximate Jacobian using finite differences; for GMM we can also pass the analytical Jacobian to gradient-based solvers by setting this 'on'
+ )
+ , silent_optimizer % run minimization of moments distance silently without displaying results or saving files in between
+
+ % Numerical algorithms options
% , aim_solver % Use AIM algorithm to compute perturbation approximation
+ % , k_order_solver % use k_order_solver in higher order perturbation approximations
% , dr=default % method used to compute the decision rule; possible values are DEFAULT, CYCLE_REDUCTION, LOGARITHMIC_REDUCTION
% , dr_cycle_reduction_tol = 1e-7 % convergence criterion used in the cycle reduction algorithm
+ % , dr_logarithmic_reduction_tol = 1e-12 % convergence criterion used in the logarithmic reduction algorithm
% , dr_logarithmic_reduction_maxiter = 100 % maximum number of iterations used in the logarithmic reduction algorithm
- % , dr_logarithmic_reduction_tol = 1e-12 % convergence criterion used in the cycle reduction algorithm
- % , k_order_solver % use k_order_solver in higher order perturbation approximations
% , lyapunov = DEFAULT % algorithm used to solve lyapunov equations; possible values are DEFAULT, FIXED_POINT, DOUBLING, SQUARE_ROOT_SOLVER
% , lyapunov_complex_threshold = 1e-15 % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
% , lyapunov_fixed_point_tol = 1e-10 % convergence criterion used in the fixed point Lyapunov solver
% , lyapunov_doubling_tol = 1e-16 % convergence criterion used in the doubling algorithm
% , sylvester = default % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
% , sylvester_fixed_point_tol = 1e-12 % convergence criterion used in the fixed point Sylvester solver
- % , qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems [IS THIS CORRET @wmutschl]
+ % , qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
% , qz_zero_threshold = 1e-6 % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
- @#if mommethod == "GMM"
- , analytic_standard_errors
- @#endif
+ % , schur_vec_tol=1e-11 % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix
+ % , mode_check % plot the target function for values around the computed minimum for each estimated parameter in turn
+ % , mode_check_neighbourhood_size = 5 % width of the window (expressed in percentage deviation) around the computed minimum to be displayed on the diagnostic plots
+ % , mode_check_symmetric_plots=1 % ensure that the check plots are symmetric around the minimum
+ % , mode_check_number_of_points = 20 % number of points around the minimum where the target function is evaluated (for each parameter)
);
@#endfor
diff --git a/tests/estimation/method_of_moments/RBC_Andreasen_Data_2.mat b/tests/estimation/method_of_moments/RBC/RBC_Andreasen_Data_2.mat
similarity index 100%
rename from tests/estimation/method_of_moments/RBC_Andreasen_Data_2.mat
rename to tests/estimation/method_of_moments/RBC/RBC_Andreasen_Data_2.mat
diff --git a/tests/estimation/method_of_moments/RBC/RBC_MoM_Andreasen.mod b/tests/estimation/method_of_moments/RBC/RBC_MoM_Andreasen.mod
new file mode 100644
index 000000000..b9191a1f7
--- /dev/null
+++ b/tests/estimation/method_of_moments/RBC/RBC_MoM_Andreasen.mod
@@ -0,0 +1,230 @@
+% Tests SMM and GMM routines
+%
+% 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 .
+% =========================================================================
+
+% Define testscenario
+@#define orderApp = 2
+@#define estimParams = 1
+
+% Note that we will set the numerical optimization tolerance levels very large to speed up the testsuite
+@#define optimizer = 13
+
+
+@#include "RBC_MoM_common.inc"
+
+shocks;
+var u_a; stderr 0.0072;
+end;
+
+varobs c iv n;
+
+
+@#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
+
+% Simulate data
+%stoch_simul(order=@{orderApp},pruning,nodisplay,nomoments,periods=500);
+%save('RBC_MoM_data_@{orderApp}.mat', options_.varobs{:} );
+%pause(1);
+
+
+estimated_params_init(use_calibration);
+end;
+
+
+%--------------------------------------------------------------------------
+% 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);
+
+c*c(-3);
+n*n(-3);
+iv*iv(-3);
+
+c*c(-5);
+n*n(-5);
+iv*iv(-5);
+end;
+
+% get indices in declaration order
+ic = strmatch('c', M_.endo_names,'exact');
+iiv = strmatch('iv', M_.endo_names,'exact');
+in = strmatch('n', M_.endo_names,'exact');
+% first entry: number of variable in declaration order
+% second entry: lag
+% third entry: power
+
+matched_moments_ = {
+ [ic ] [0 ], [1 ];
+ [in ] [0 ], [1 ];
+ [iiv ] [0 ], [1 ];
+
+ [ic ic ] [0 0], [1 1];
+ [ic iiv] [0 0], [1 1];
+ %[ic in ] [0 0], [1 1];
+ %[iiv ic ] [0 0], [1 1];
+ [in iiv] [0 0], [1 1];
+ [iiv iiv] [0 0], [1 1];
+ [ic in] [0 0], [1 1];
+ %[in iiv] [0 0], [1 1];
+ [in in ] [0 0], [1 1];
+
+ [ic ic ] [0 -1], [1 1];
+ [in in ] [0 -1], [1 1];
+ [iiv iiv] [0 -1], [1 1];
+
+ [ic ic ] [0 -3], [1 1];
+ [in in ] [0 -3], [1 1];
+ [iiv iiv] [0 -3], [1 1];
+
+ [ic ic ] [0 -5], [1 1];
+ [in in ] [0 -5], [1 1];
+ [iiv iiv] [0 -5], [1 1];
+
+};
+
+if ~isequal(M_.matched_moments,matched_moments_)
+ error('Translation to matched_moments-block failed')
+end
+
+
+method_of_moments(
+ % Necessery options
+ mom_method = GMM % method of moments method; possible values: GMM|SMM
+ , datafile = 'RBC_Andreasen_Data_2.mat' % name of filename with data
+
+ % Options for both GMM and SMM
+ % , bartlett_kernel_lag = 20 % bandwith in optimal weighting matrix
+ , order = @{orderApp} % order of Taylor approximation in perturbation
+ % , penalized_estimator % include deviation from prior mean as additional moment restriction and use prior precision as weight
+ % , pruning % use pruned state space system at higher-order
+ % , verbose % display and store intermediate estimation results
+ , weighting_matrix = ['DIAGONAL','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']
+ % , weighting_matrix_scaling_factor=1 % scaling of weighting matrix in objective function
+ , se_tolx=1e-6 % step size for numerical computation of standard errors
+
+ % Options for SMM
+ % , burnin=500 % number of periods dropped at beginning of simulation
+ % , bounded_shock_support % trim shocks in simulation to +- 2 stdev
+ % , seed = 24051986 % seed used in simulations
+ % , simulation_multiple = 5 % multiple of the data length used for simulation
+
+ % Options for GMM
+ % , analytic_standard_errors % compute standard errors using analytical derivatives
+
+ % General options
+ % , dirname = 'MM' % directory in which to store estimation output
+ % , graph_format = EPS % specify the file format(s) for graphs saved to disk
+ % , nodisplay % do not display the graphs, but still save them to disk
+ % , nograph % do not create graphs (which implies that they are not saved to the disk nor displayed)
+ % , noprint % do not print stuff to console
+ % , plot_priors = 1 % control plotting of priors
+ % , prior_trunc = 1e-10 % probability of extreme values of the prior density that is ignored when computing bounds for the parameters
+ , TeX % print TeX tables and graphics
+
+ % Data and model options
+ % , first_obs = 501 % number of first observation
+ % , logdata % if data is already in logs
+ % , nobs = 250 % number of observations
+ % , prefilter=0 % demean each data series by its empirical mean and use centered moments
+ % , xls_sheet = data % name/number of sheet with data in Excel
+ % , xls_range = B2:D200 % range of data in Excel sheet
+
+ % Optimization options that can be set by the user in the mod file, otherwise default values are provided
+ % , huge_number=1e7 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+ , mode_compute = 3 % specifies the optimizer for minimization of moments distance
+ , additional_optimizer_steps = [13] % vector of additional mode-finders run after mode_compute
+ % optim: a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute, some exemplary common options:
+ , optim = ('TolFun' , 1D-6 % termination tolerance on the function value, a positive scalar
+ ,'TolX' , 1e-6 % termination tolerance on x, a positive scalar
+ % ,'MaxIter' , 3000 % maximum number of iterations allowed, a positive integer
+ % ,'MaxFunEvals' , 1D6 % maximum number of function evaluations allowed, a positive integer
+ % ,'UseParallel' , 1 % when true (and supported by optimizer) solver estimates gradients in parallel (using Matlab/Octave's parallel toolbox)
+ % ,'Jacobian' , 'off' % when 'off' gradient-based solvers approximate Jacobian using finite differences; for GMM we can also pass the analytical Jacobian to gradient-based solvers by setting this 'on'
+ )
+ % , silent_optimizer % run minimization of moments distance silently without displaying results or saving files in between
+
+ % Numerical algorithms options
+ % , aim_solver % Use AIM algorithm to compute perturbation approximation
+ % , k_order_solver % use k_order_solver in higher order perturbation approximations
+ % , dr=default % method used to compute the decision rule; possible values are DEFAULT, CYCLE_REDUCTION, LOGARITHMIC_REDUCTION
+ % , dr_cycle_reduction_tol = 1e-7 % convergence criterion used in the cycle reduction algorithm
+ % , dr_logarithmic_reduction_tol = 1e-12 % convergence criterion used in the logarithmic reduction algorithm
+ % , dr_logarithmic_reduction_maxiter = 100 % maximum number of iterations used in the logarithmic reduction algorithm
+ % , lyapunov = DEFAULT % algorithm used to solve lyapunov equations; possible values are DEFAULT, FIXED_POINT, DOUBLING, SQUARE_ROOT_SOLVER
+ % , lyapunov_complex_threshold = 1e-15 % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
+ % , lyapunov_fixed_point_tol = 1e-10 % convergence criterion used in the fixed point Lyapunov solver
+ % , lyapunov_doubling_tol = 1e-16 % convergence criterion used in the doubling algorithm
+ % , sylvester = default % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
+ % , sylvester_fixed_point_tol = 1e-12 % convergence criterion used in the fixed point Sylvester solver
+ % , qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
+ % , qz_zero_threshold = 1e-6 % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
+ % , schur_vec_tol=1e-11 % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix
+ , mode_check % plot the target function for values around the computed minimum for each estimated parameter in turn
+ % , mode_check_neighbourhood_size = 5 % width of the window (expressed in percentage deviation) around the computed minimum to be displayed on the diagnostic plots
+ % , mode_check_symmetric_plots=1 % ensure that the check plots are symmetric around the minimum
+ % , mode_check_number_of_points = 20 % number of points around the minimum where the target function is evaluated (for each parameter)
+ );
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/RBC_MoM_SMM_ME.mod b/tests/estimation/method_of_moments/RBC/RBC_MoM_SMM_ME.mod
similarity index 65%
rename from tests/estimation/method_of_moments/RBC_MoM_SMM_ME.mod
rename to tests/estimation/method_of_moments/RBC/RBC_MoM_SMM_ME.mod
index 964862670..b407cd28a 100644
--- a/tests/estimation/method_of_moments/RBC_MoM_SMM_ME.mod
+++ b/tests/estimation/method_of_moments/RBC/RBC_MoM_SMM_ME.mod
@@ -1,3 +1,5 @@
+% =========================================================================
+% Copyright (C) 2020-2021 Dynare Team
%
% This file is part of Dynare.
%
@@ -137,30 +139,31 @@ end
@#for mommethod in ["SMM"]
method_of_moments(
- % Necessery options
- mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM
- , datafile = 'RBC_MoM_data_@{orderApp}.mat' % name of filename with data
+ % Necessery options
+ mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM
+ , datafile = 'RBC_MoM_data_@{orderApp}.mat' % name of filename with data
- % Options for both GMM and SMM
+ % Options for both GMM and SMM
% , bartlett_kernel_lag = 20 % bandwith in optimal weighting matrix
, order = @{orderApp} % order of Taylor approximation in perturbation
- % , penalized_estimator % use penalized optimization
+ % , penalized_estimator % include deviation from prior mean as additional moment restriction and use prior precision as weight
, pruning % use pruned state space system at higher-order
% , verbose % display and store intermediate estimation results
- , weighting_matrix = ['identity_matrix'] % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
- , weighting_matrix_scaling_factor = 10
- , burnin=250
- %, additional_optimizer_steps = [4] % vector of additional mode-finders run after mode_compute
- % , prefilter=0 % demean each data series by its empirical mean and use centered moments
- %
- % Options for SMM
+ , weighting_matrix = ['identity_matrix'] % 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']
+ , weighting_matrix_scaling_factor=10 % scaling of weighting matrix in objective function
+ % , se_tolx=1e-6 % step size for numerical computation of standard errors
+
+ % Options for SMM
+ , burnin=250 % number of periods dropped at beginning of simulation
% , bounded_shock_support % trim shocks in simulation to +- 2 stdev
- % , drop = 500 % number of periods dropped at beginning of simulation
% , seed = 24051986 % seed used in simulations
% , simulation_multiple = 5 % multiple of the data length used for simulation
- %
- % General options
- %, dirname = 'MM' % directory in which to store estimation output
+
+ % Options for GMM
+ % , analytic_standard_errors % compute standard errors using analytical derivatives
+
+ % General options
+ % , dirname = 'MM' % directory in which to store estimation output
% , graph_format = EPS % specify the file format(s) for graphs saved to disk
% , nodisplay % do not display the graphs, but still save them to disk
% , nograph % do not create graphs (which implies that they are not saved to the disk nor displayed)
@@ -168,41 +171,49 @@ end
% , plot_priors = 1 % control plotting of priors
% , prior_trunc = 1e-10 % probability of extreme values of the prior density that is ignored when computing bounds for the parameters
% , TeX % print TeX tables and graphics
- %
- % Data and model options
- %, first_obs = 501 % number of first observation
- % , logdata % if loglinear is set, this option is necessary if the user provides data already in logs, otherwise the log transformation will be applied twice (this may result in complex data)
- % , loglinear % computes a log-linear approximation of the model instead of a linear approximation
- %, nobs = 500 % number of observations
- % , xls_sheet = willi % name of sheet with data in Excel
+
+ % Data and model options
+ % , first_obs = 501 % number of first observation
+ % , logdata % if data is already in logs
+ % , nobs = 250 % number of observations
+ % , prefilter=0 % demean each data series by its empirical mean and use centered moments
+ % , xls_sheet = data % name/number of sheet with data in Excel
% , xls_range = B2:D200 % range of data in Excel sheet
- %
- % Optimization options that can be set by the user in the mod file, otherwise default values are provided
- % , analytic_derivation % uses analytic derivatives to compute standard errors for GMM
- %, huge_number=1D10 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
- , mode_compute = @{optimizer} % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
- %, optim = ('TolFun', 1e-3
- % ,'TolX', 1e-5
- % ) % 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
- % , tolf = 1e-5 % convergence criterion on function value for numerical differentiation
- % , tolx = 1e-6 % convergence criterion on funciton input for numerical differentiation
- %
- % % Numerical algorithms options
+
+ % Optimization options that can be set by the user in the mod file, otherwise default values are provided
+ % , huge_number=1e7 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+ , mode_compute = @{optimizer} % specifies the optimizer for minimization of moments distance
+ %, additional_optimizer_steps = [1 2 3 4] % vector of additional mode-finders run after mode_compute
+ % optim: a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute, some exemplary common options:
+ % , optim = ('TolFun' , 1e-6 % termination tolerance on the function value, a positive scalar
+ % ,'TolX' , 1e-6 % termination tolerance on x, a positive scalar
+ % ,'MaxIter' , 3000 % maximum number of iterations allowed, a positive integer
+ % ,'MaxFunEvals' , 1D6 % maximum number of function evaluations allowed, a positive integer
+ % ,'UseParallel' , 1 % when true (and supported by optimizer) solver estimates gradients in parallel (using Matlab/Octave's parallel toolbox)
+ % ,'Jacobian' , 'off' % when 'off' gradient-based solvers approximate Jacobian using finite differences; for GMM we can also pass the analytical Jacobian to gradient-based solvers by setting this 'on'
+ % )
+ % , silent_optimizer % run minimization of moments distance silently without displaying results or saving files in between
+
+ % Numerical algorithms options
% , aim_solver % Use AIM algorithm to compute perturbation approximation
+ % , k_order_solver % use k_order_solver in higher order perturbation approximations
% , dr=default % method used to compute the decision rule; possible values are DEFAULT, CYCLE_REDUCTION, LOGARITHMIC_REDUCTION
% , dr_cycle_reduction_tol = 1e-7 % convergence criterion used in the cycle reduction algorithm
+ % , dr_logarithmic_reduction_tol = 1e-12 % convergence criterion used in the logarithmic reduction algorithm
% , dr_logarithmic_reduction_maxiter = 100 % maximum number of iterations used in the logarithmic reduction algorithm
- % , dr_logarithmic_reduction_tol = 1e-12 % convergence criterion used in the cycle reduction algorithm
- % , k_order_solver % use k_order_solver in higher order perturbation approximations
% , lyapunov = DEFAULT % algorithm used to solve lyapunov equations; possible values are DEFAULT, FIXED_POINT, DOUBLING, SQUARE_ROOT_SOLVER
% , lyapunov_complex_threshold = 1e-15 % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
% , lyapunov_fixed_point_tol = 1e-10 % convergence criterion used in the fixed point Lyapunov solver
% , lyapunov_doubling_tol = 1e-16 % convergence criterion used in the doubling algorithm
% , sylvester = default % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
% , sylvester_fixed_point_tol = 1e-12 % convergence criterion used in the fixed point Sylvester solver
- % , qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems [IS THIS CORRET @wmutschl]
+ % , qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
% , qz_zero_threshold = 1e-6 % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
+ % , schur_vec_tol=1e-11 % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix
+ % , mode_check % plot the target function for values around the computed minimum for each estimated parameter in turn
+ % , mode_check_neighbourhood_size = 5 % width of the window (expressed in percentage deviation) around the computed minimum to be displayed on the diagnostic plots
+ % , mode_check_symmetric_plots=1 % ensure that the check plots are symmetric around the minimum
+ % , mode_check_number_of_points = 20 % number of points around the minimum where the target function is evaluated (for each parameter)
);
@#endfor
diff --git a/tests/estimation/method_of_moments/RBC_MoM_common.inc b/tests/estimation/method_of_moments/RBC/RBC_MoM_common.inc
similarity index 71%
rename from tests/estimation/method_of_moments/RBC_MoM_common.inc
rename to tests/estimation/method_of_moments/RBC/RBC_MoM_common.inc
index d480e35c0..330dd2fcf 100644
--- a/tests/estimation/method_of_moments/RBC_MoM_common.inc
+++ b/tests/estimation/method_of_moments/RBC/RBC_MoM_common.inc
@@ -2,7 +2,23 @@
% Andreasen, Fernández-Villaverde, Rubio-RamÃrez (2018): "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications", Review of Economic Studies, 85(1):1-49.
% Adapted by Willi Mutschler (@wmutschl, willi@mutschler.eu)
% =========================================================================
-% Copyright (C) 2020 Dynare Team
+% 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 .
+% =========================================================================
var k $K$
c $C$
diff --git a/tests/estimation/method_of_moments/RBC/RBC_MoM_optimizer.mod b/tests/estimation/method_of_moments/RBC/RBC_MoM_optimizer.mod
new file mode 100644
index 000000000..08b1ed882
--- /dev/null
+++ b/tests/estimation/method_of_moments/RBC/RBC_MoM_optimizer.mod
@@ -0,0 +1,146 @@
+% Test optimizers
+%
+% 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 .
+% =========================================================================
+% TO DO
+% [ ] fix optimizers 11 and 12;
+% note that 12 and 102 require GADS_Toolbox which is not available on servers, but need to be tested locally
+
+% Define testscenario
+@#define orderApp = 2
+
+% Note that we will set the numerical optimization tolerance levels very large to speed up the testsuite
+
+@#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;
+
+% reduce options to speed up testsuite
+options_.newrat.maxiter = 10;
+options_.newrat.tolerance.f = 1e-2;
+options_.newrat.tolerance.f_analytic = 1e-2;
+
+options_.mh_jscale = 0.6;
+options_.gmhmaxlik.iterations=1;
+options_.gmhmaxlik.number=2000;
+options_.gmhmaxlik.nclimb=2000;
+options_.gmhmaxlik.nscale=2000;
+options_.gmhmaxlik.target=0.5;
+
+options_.solveopt.MaxIter=300;
+options_.solveopt.LBGradientStep=1e-3;
+options_.solveopt.TolFun = 1e-3;
+options_.solveopt.TolX = 1e-3;
+options_.solveopt.TolXConstraint=1e-3;
+
+
+@#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;
+ @#define OPTIMIZERS = [1, 2, 3, 4, 5, 7, 8, 9, 10, 13, 101]
+ @#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;
+ @#define OPTIMIZERS = [1, 2, 3, 4, 7, 8, 9, 10, 13, 101]
+ @#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;
+ @#define OPTIMIZERS = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 13, 101]
+ @#endif
+
+ estimated_params_init(use_calibration);
+ end;
+
+ @#for optimizer in OPTIMIZERS
+ 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 = @{orderApp} % 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']
+ , nograph % do not create graphs (which implies that they are not saved to the disk nor displayed)
+ , mode_compute = @{optimizer} % specifies the optimizer for minimization of moments distance
+ @#if optimizer == 102
+ , optim = ('TolFun' , 1D-3 % termination tolerance on the function value, a positive scalar
+ ,'MaxIter' , 300 % maximum number of iterations allowed, a positive integer
+ ,'MaxFunEvals' , 1D3 % maximum number of function evaluations allowed, a positive integer
+ )
+ @#else
+ , optim = ('TolFun' , 1D-3 % termination tolerance on the function value, a positive scalar
+ ,'TolX' , 1e-3 % termination tolerance on x, a positive scalar
+ ,'MaxIter' , 300 % maximum number of iterations allowed, a positive integer
+ ,'MaxFunEvals' , 1D3 % maximum number of function evaluations allowed, a positive integer
+ )
+ @#endif
+ %, silent_optimizer % run minimization of moments distance silently without displaying results or saving files in between
+ );
+ @#endfor
+@#endfor
+
+
+
diff --git a/tests/estimation/method_of_moments/RBC_MoM_prefilter.mod b/tests/estimation/method_of_moments/RBC/RBC_MoM_prefilter.mod
similarity index 61%
rename from tests/estimation/method_of_moments/RBC_MoM_prefilter.mod
rename to tests/estimation/method_of_moments/RBC/RBC_MoM_prefilter.mod
index 7fb29f8ab..326badcd5 100644
--- a/tests/estimation/method_of_moments/RBC_MoM_prefilter.mod
+++ b/tests/estimation/method_of_moments/RBC/RBC_MoM_prefilter.mod
@@ -1,6 +1,6 @@
% Tests SMM and GMM routines with prefilter, explicit initialization, and estimated_params_init(use_calibration);
%
-% Copyright (C) 2020 Dynare Team
+% Copyright (C) 2020-2021 Dynare Team
%
% This file is part of Dynare.
%
@@ -110,31 +110,31 @@ save('test_matrix.mat','weighting_matrix')
@#for mommethod in ["GMM", "SMM"]
method_of_moments(
- % Necessery options
- mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM
- , datafile = 'RBC_MoM_data_@{orderApp}.mat' % name of filename with data
+ % Necessery options
+ mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM
+ , datafile = 'RBC_MoM_data_@{orderApp}.mat' % name of filename with data
- % Options for both GMM and SMM
+ % Options for both GMM and SMM
% , bartlett_kernel_lag = 20 % bandwith in optimal weighting matrix
, order = @{orderApp} % order of Taylor approximation in perturbation
- % , penalized_estimator % use penalized optimization
+ % , penalized_estimator % include deviation from prior mean as additional moment restriction and use prior precision as weight
, pruning % use pruned state space system at higher-order
% , verbose % display and store intermediate estimation results
-% , weighting_matrix = 'test_matrix.mat' % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
- , weighting_matrix =['test_matrix.mat','optimal']
- %, weighting_matrix = optimal % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
- %, additional_optimizer_steps = [4] % vector of additional mode-finders run after mode_compute
- , prefilter=1 % demean each data series by its empirical mean and use centered moments
- , se_tolx=1e-5
- %
- % Options for SMM
+ , weighting_matrix = ['test_matrix.mat','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']
+ %, weighting_matrix_scaling_factor=1 % scaling of weighting matrix in objective function
+ , se_tolx=1e-5 % step size for numerical computation of standard errors
+
+ % Options for SMM
+ , burnin=500 % number of periods dropped at beginning of simulation
% , bounded_shock_support % trim shocks in simulation to +- 2 stdev
- , burnin = 500 % number of periods dropped at beginning of simulation
% , seed = 24051986 % seed used in simulations
% , simulation_multiple = 5 % multiple of the data length used for simulation
- %
- % General options
- %, dirname = 'MM' % directory in which to store estimation output
+
+ % Options for GMM
+ % , analytic_standard_errors % compute standard errors using analytical derivatives
+
+ % General options
+ % , dirname = 'MM' % directory in which to store estimation output
% , graph_format = EPS % specify the file format(s) for graphs saved to disk
% , nodisplay % do not display the graphs, but still save them to disk
% , nograph % do not create graphs (which implies that they are not saved to the disk nor displayed)
@@ -142,38 +142,49 @@ save('test_matrix.mat','weighting_matrix')
% , plot_priors = 1 % control plotting of priors
% , prior_trunc = 1e-10 % probability of extreme values of the prior density that is ignored when computing bounds for the parameters
% , TeX % print TeX tables and graphics
- %
- % Data and model options
- %, first_obs = 501 % number of first observation
- % , logdata % if loglinear is set, this option is necessary if the user provides data already in logs, otherwise the log transformation will be applied twice (this may result in complex data)
- % , loglinear % computes a log-linear approximation of the model instead of a linear approximation
- %, nobs = 500 % number of observations
- % , xls_sheet = willi % name of sheet with data in Excel
+
+ % Data and model options
+ % , first_obs = 501 % number of first observation
+ % , logdata % if data is already in logs
+ , nobs = 250 % number of observations
+ , prefilter=1 % demean each data series by its empirical mean and use centered moments
+
+ % , xls_sheet = data % name/number of sheet with data in Excel
% , xls_range = B2:D200 % range of data in Excel sheet
- %
- % Optimization options that can be set by the user in the mod file, otherwise default values are provided
- % , analytic_derivation % uses analytic derivatives to compute standard errors for GMM
- %, huge_number=1D10 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
- , mode_compute = @{optimizer} % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
- %, optim = ('TolFun', 1e-3
- % ,'TolX', 1e-5
- % ) % 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
- %
- % % Numerical algorithms options
+
+ % Optimization options that can be set by the user in the mod file, otherwise default values are provided
+ % , huge_number=1e7 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
+ , mode_compute = @{optimizer} % specifies the optimizer for minimization of moments distance
+ %, additional_optimizer_steps = [7] % vector of additional mode-finders run after mode_compute
+ % optim: a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute, some exemplary common options:
+ % , optim = ('TolFun' , 1e-6 % termination tolerance on the function value, a positive scalar
+ % ,'TolX' , 1e-6 % termination tolerance on x, a positive scalar
+ % ,'MaxIter' , 3000 % maximum number of iterations allowed, a positive integer
+ % ,'MaxFunEvals' , 1D6 % maximum number of function evaluations allowed, a positive integer
+ % ,'UseParallel' , 1 % when true (and supported by optimizer) solver estimates gradients in parallel (using Matlab/Octave's parallel toolbox)
+ % ,'Jacobian' , 'off' % when 'off' gradient-based solvers approximate Jacobian using finite differences; for GMM we can also pass the analytical Jacobian to gradient-based solvers by setting this 'on'
+ % )
+ % , silent_optimizer % run minimization of moments distance silently without displaying results or saving files in between
+
+ % Numerical algorithms options
% , aim_solver % Use AIM algorithm to compute perturbation approximation
+ % , k_order_solver % use k_order_solver in higher order perturbation approximations
% , dr=default % method used to compute the decision rule; possible values are DEFAULT, CYCLE_REDUCTION, LOGARITHMIC_REDUCTION
% , dr_cycle_reduction_tol = 1e-7 % convergence criterion used in the cycle reduction algorithm
+ % , dr_logarithmic_reduction_tol = 1e-12 % convergence criterion used in the logarithmic reduction algorithm
% , dr_logarithmic_reduction_maxiter = 100 % maximum number of iterations used in the logarithmic reduction algorithm
- % , dr_logarithmic_reduction_tol = 1e-12 % convergence criterion used in the cycle reduction algorithm
- % , k_order_solver % use k_order_solver in higher order perturbation approximations
% , lyapunov = DEFAULT % algorithm used to solve lyapunov equations; possible values are DEFAULT, FIXED_POINT, DOUBLING, SQUARE_ROOT_SOLVER
% , lyapunov_complex_threshold = 1e-15 % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
% , lyapunov_fixed_point_tol = 1e-10 % convergence criterion used in the fixed point Lyapunov solver
% , lyapunov_doubling_tol = 1e-16 % convergence criterion used in the doubling algorithm
% , sylvester = default % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
% , sylvester_fixed_point_tol = 1e-12 % convergence criterion used in the fixed point Sylvester solver
- % , qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems [IS THIS CORRET @wmutschl]
+ % , qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
% , qz_zero_threshold = 1e-6 % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
+ % , schur_vec_tol=1e-11 % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix
+ % , mode_check % plot the target function for values around the computed minimum for each estimated parameter in turn
+ % , mode_check_neighbourhood_size = 5 % width of the window (expressed in percentage deviation) around the computed minimum to be displayed on the diagnostic plots
+ % , mode_check_symmetric_plots=1 % ensure that the check plots are symmetric around the minimum
+ % , mode_check_number_of_points = 20 % number of points around the minimum where the target function is evaluated (for each parameter)
);
@#endfor
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/RBC/RBC_MoM_steady_helper.m b/tests/estimation/method_of_moments/RBC/RBC_MoM_steady_helper.m
new file mode 100644
index 000000000..9c43619d7
--- /dev/null
+++ b/tests/estimation/method_of_moments/RBC/RBC_MoM_steady_helper.m
@@ -0,0 +1,40 @@
+% =========================================================================
+% 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 .
+% =========================================================================
+function [N, info]= RBC_MoM_steady_helper(THETA,ETAl,ETAc,BETTA,B,C_O_N,W)
+info=0;
+if ~isreal(C_O_N)
+ info=1;
+ N=NaN;
+ return;
+end
+if ETAc == 1 && ETAl == 1
+ N = (1-BETTA*B)*(C_O_N*(1-B))^-1*W/THETA/(1+(1-BETTA*B)*(C_O_N*(1-B))^-1*W/THETA);
+else
+ % No closed-form solution use a fixed-point algorithm
+ N0 = 1/3;
+ try
+ [N, ~, exitflag] = fsolve(@(N) THETA*(1-N)^(-ETAl)*N^ETAc - (1-BETTA*B)*(C_O_N*(1-B))^(-ETAc)*W, N0,optimset('Display','off','TolX',1e-12,'TolFun',1e-12));
+ if exitflag<1
+ info=1;
+ end
+ catch
+ N=NaN;
+ info=1;
+ end
+end
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/RBC_MoM_Andreasen.mod b/tests/estimation/method_of_moments/RBC_MoM_Andreasen.mod
deleted file mode 100644
index ae6984005..000000000
--- a/tests/estimation/method_of_moments/RBC_MoM_Andreasen.mod
+++ /dev/null
@@ -1,227 +0,0 @@
-% Tests SMM and GMM routines
-%
-% Copyright (C) 2020 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 .
-% =========================================================================
-
-% Define testscenario
-@#define orderApp = 2
-@#define estimParams = 1
-
-% Note that we will set the numerical optimization tolerance levels very large to speed up the testsuite
-@#define optimizer = 13
-
-
-@#include "RBC_MoM_common.inc"
-
-shocks;
-var u_a; stderr 0.0072;
-end;
-
-varobs c iv n;
-
-
-@#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
-
-% Simulate data
-%stoch_simul(order=@{orderApp},pruning,nodisplay,nomoments,periods=500);
-%save('RBC_MoM_data_@{orderApp}.mat', options_.varobs{:} );
-%pause(1);
-
-
-estimated_params_init(use_calibration);
-end;
-
-
-%--------------------------------------------------------------------------
-% 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);
-
-c*c(-3);
-n*n(-3);
-iv*iv(-3);
-
-c*c(-5);
-n*n(-5);
-iv*iv(-5);
-end;
-
-% get indices in declaration order
-ic = strmatch('c', M_.endo_names,'exact');
-iiv = strmatch('iv', M_.endo_names,'exact');
-in = strmatch('n', M_.endo_names,'exact');
-% first entry: number of variable in declaration order
-% second entry: lag
-% third entry: power
-
-matched_moments_ = {
- [ic ] [0 ], [1 ];
- [in ] [0 ], [1 ];
- [iiv ] [0 ], [1 ];
-
- [ic ic ] [0 0], [1 1];
- [ic iiv] [0 0], [1 1];
- %[ic in ] [0 0], [1 1];
- %[iiv ic ] [0 0], [1 1];
- [in iiv] [0 0], [1 1];
- [iiv iiv] [0 0], [1 1];
- [ic in] [0 0], [1 1];
- %[in iiv] [0 0], [1 1];
- [in in ] [0 0], [1 1];
-
- [ic ic ] [0 -1], [1 1];
- [in in ] [0 -1], [1 1];
- [iiv iiv] [0 -1], [1 1];
-
- [ic ic ] [0 -3], [1 1];
- [in in ] [0 -3], [1 1];
- [iiv iiv] [0 -3], [1 1];
-
- [ic ic ] [0 -5], [1 1];
- [in in ] [0 -5], [1 1];
- [iiv iiv] [0 -5], [1 1];
-
-};
-
-if ~isequal(M_.matched_moments,matched_moments_)
- error('Translation to matched_moments-block failed')
-end
-
-
- method_of_moments(
- % Necessery options
- mom_method = GMM % method of moments method; possible values: GMM|SMM
- , datafile = 'RBC_Andreasen_Data_2.mat' % name of filename with data
-
- % Options for both GMM and SMM
- %, bartlett_kernel_lag = 20 % bandwith in optimal weighting matrix
- , order = @{orderApp} % order of Taylor approximation in perturbation
- %, penalized_estimator % use penalized optimization
- %, pruning % use pruned state space system at higher-order
- %, verbose % display and store intermediate estimation results
- , weighting_matrix = ['DIAGONAL','OPTIMAL'] % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
- %, weighting_matrix_scaling_factor=1
- , additional_optimizer_steps = [13] % vector of additional mode-finders run after mode_compute
- %, prefilter=0 % demean each data series by its empirical mean and use centered moments
- %
- % Options for SMM
- %, bounded_shock_support % trim shocks in simulation to +- 2 stdev
- %, drop = 500 % number of periods dropped at beginning of simulation
- %, seed = 24051986 % seed used in simulations
- %, simulation_multiple = 5 % multiple of the data length used for simulation
- %, burnin = 200
- %
- % General options
- %, dirname = 'MM' % directory in which to store estimation output
- %, graph_format = EPS % specify the file format(s) for graphs saved to disk
- %, nodisplay % do not display the graphs, but still save them to disk
- %, nograph % do not create graphs (which implies that they are not saved to the disk nor displayed)
- %, noprint % do not print stuff to console
- %, plot_priors = 1 % control plotting of priors
- %, prior_trunc = 1e-10 % probability of extreme values of the prior density that is ignored when computing bounds for the parameters
- , TeX % print TeX tables and graphics
- %
- % Data and model options
- %, first_obs = 501 % number of first observation
- %, logdata % if loglinear is set, this option is necessary if the user provides data already in logs, otherwise the log transformation will be applied twice (this may result in complex data)
- %, loglinear % computes a log-linear approximation of the model instead of a linear approximation
- %, nobs = 50 % number of observations
- % , xls_sheet = willi % name of sheet with data in Excel
- % , xls_range = B2:D200 % range of data in Excel sheet
- %
- % Optimization options that can be set by the user in the mod file, otherwise default values are provided
- %, analytic_derivation % uses analytic derivatives to compute standard errors for GMM
- %, huge_number=1D10 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
- , mode_compute = 13 % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
- , optim = ('TolFun', 1D-6
- ,'TolX', 1D-6
- ) % 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
- , se_tolx = 1e-6 % convergence criterion on funciton input for numerical differentiation
- %
- % % Numerical algorithms options
- %, aim_solver % Use AIM algorithm to compute perturbation approximation
- %, dr=DEFAULT % method used to compute the decision rule; possible values are DEFAULT, CYCLE_REDUCTION, LOGARITHMIC_REDUCTION
- %, dr_cycle_reduction_tol = 1e-7 % convergence criterion used in the cycle reduction algorithm
- %, dr_logarithmic_reduction_maxiter = 100 % maximum number of iterations used in the logarithmic reduction algorithm
- %, dr_logarithmic_reduction_tol = 1e-12 % convergence criterion used in the cycle reduction algorithm
- %, k_order_solver % use k_order_solver in higher order perturbation approximations
- %, lyapunov = DEFAULT % algorithm used to solve lyapunov equations; possible values are DEFAULT, FIXED_POINT, DOUBLING, SQUARE_ROOT_SOLVER
- %, lyapunov_complex_threshold = 1e-15 % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
- %, lyapunov_fixed_point_tol = 1e-10 % convergence criterion used in the fixed point Lyapunov solver
- %, lyapunov_doubling_tol = 1e-16 % convergence criterion used in the doubling algorithm
- %, sylvester = default % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
- %, sylvester_fixed_point_tol = 1e-12 % convergence criterion used in the fixed point Sylvester solver
- %, qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems [IS THIS CORRET @wmutschl]
- %, qz_zero_threshold = 1e-6 % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
- , mode_check
- %, mode_check_neighbourhood_size=0.5
- %, mode_check_symmetric_plots=0
- %, mode_check_number_of_points=25
- );
-
-
-
diff --git a/tests/estimation/method_of_moments/RBC_MoM_steady_helper.m b/tests/estimation/method_of_moments/RBC_MoM_steady_helper.m
deleted file mode 100644
index b495e27a1..000000000
--- a/tests/estimation/method_of_moments/RBC_MoM_steady_helper.m
+++ /dev/null
@@ -1,22 +0,0 @@
-function [N, info]= RBC_MoM_steady_helper(THETA,ETAl,ETAc,BETTA,B,C_O_N,W)
-info=0;
-if ~isreal(C_O_N)
- info=1;
- N=NaN;
- return;
-end
-if ETAc == 1 && ETAl == 1
- N = (1-BETTA*B)*(C_O_N*(1-B))^-1*W/THETA/(1+(1-BETTA*B)*(C_O_N*(1-B))^-1*W/THETA);
-else
- % No closed-form solution use a fixed-point algorithm
- N0 = 1/3;
- try
- [N, ~, exitflag] = fsolve(@(N) THETA*(1-N)^(-ETAl)*N^ETAc - (1-BETTA*B)*(C_O_N*(1-B))^(-ETAc)*W, N0,optimset('Display','off','TolX',1e-12,'TolFun',1e-12));
- if exitflag<1
- info=1;
- end
- catch
- N=NaN;
- info=1;
- end
-end
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/RBC_MoM_steadystate.m b/tests/estimation/method_of_moments/RBC_MoM_steadystate.m
deleted file mode 100644
index ba4ef9240..000000000
--- a/tests/estimation/method_of_moments/RBC_MoM_steadystate.m
+++ /dev/null
@@ -1,74 +0,0 @@
-% By Willi Mutschler, September 26, 2016. Email: willi@mutschler.eu
-function [ys,params,check] = RBCmodel_steadystate(ys,exo,M_,options_)
-%% Step 0: initialize indicator and set options for numerical solver
-check = 0;
-options = optimset('Display','off','TolX',1e-12,'TolFun',1e-12);
-params = M_.params;
-
-%% Step 1: read out parameters to access them with their name
-for ii = 1:M_.param_nbr
- eval([ M_.param_names{ii} ' = M_.params(' int2str(ii) ');']);
-end
-
-%% Step 2: Check parameter restrictions
-if ETAc*ETAl<1 % parameter violates restriction (here it is artifical)
- check=1; %set failure indicator
- return; %return without updating steady states
-end
-
-%% Step 3: Enter model equations here
-A = 1;
-RK = 1/BETTA - (1-DELTA);
-K_O_N = (RK/(A*(1-ALFA)))^(-1/ALFA);
-if K_O_N <= 0
- check = 1; % set failure indicator
- return; % return without updating steady states
-end
-W = A*ALFA*(K_O_N)^(1-ALFA);
-IV_O_N = DELTA*K_O_N;
-Y_O_N = A*K_O_N^(1-ALFA);
-C_O_N = Y_O_N - IV_O_N;
-if C_O_N <= 0
- check = 1; % set failure indicator
- return; % return without updating steady states
-end
-
-% The labor level
-if ETAc == 1 && ETAl == 1
- N = (1-BETTA*B)*(C_O_N*(1-B))^-1*W/THETA/(1+(1-BETTA*B)*(C_O_N*(1-B))^-1*W/THETA);
-else
- % No closed-form solution use a fixed-point algorithm
- N0 = 1/3;
- [N,~,exitflag] = fsolve(@(N) THETA*(1-N)^(-ETAl)*N^ETAc - (1-BETTA*B)*(C_O_N*(1-B))^(-ETAc)*W, N0,options);
- if exitflag <= 0
- check = 1; % set failure indicator
- return % return without updating steady states
- end
-end
-
-C=C_O_N*N;
-Y=Y_O_N*N;
-IV=IV_O_N*N;
-K=K_O_N*N;
-LA = (C-B*C)^(-ETAc)-BETTA*B*(C-B*C)^(-ETAc);
-
-k=log(K);
-c=log(C);
-a=log(A);
-iv=log(IV);
-y=log(Y);
-la=log(LA);
-n=log(N);
-rk=log(RK);
-w=log(W);
-%% Step 4: Update parameters and variables
-params=NaN(M_.param_nbr,1);
-for iter = 1:M_.param_nbr %update parameters set in the file
- eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
-end
-
-for ii = 1:M_.orig_endo_nbr %auxiliary variables are set automatically
- eval(['ys(' int2str(ii) ') = ' M_.endo_names{ii} ';']);
-end
-
-end