MoM: Refactor AnScho test case

Instead of checking everything in one mod file, this commit separates the checks into individual mod files that test:
- whether the translation from matched_moments works
- whether the duplicate moments are found
- whether GMM and SMM both work with different estimated_params blocks.
wip

tests/.gitignore

@@ -50,6 +50,7 @@ wsOct
 !/ep/mean_preserving_spread.m
 !/ep/rbcii_steady_state.m
 !/estimation/fsdat_simul.m
+!/estimation/method_of_moments/AnScho/AnScho_MoM_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

tests/Makefile.am

@@ -62,7 +62,13 @@ MODFILES = \
 	estimation/heteroskedastic_shocks/fs2000_het_corr.mod \
 	estimation/t_proposal/fs2000_student.mod \
 	estimation/tune_mh_jscale/fs2000.mod \
-	estimation/method_of_moments/AnScho/AnScho_MoM.mod \
+	estimation/method_of_moments/AnScho/AnScho_matched_moments.mod \
+	estimation/method_of_moments/AnScho/AnScho_GMM_estimParams0.mod \
+	estimation/method_of_moments/AnScho/AnScho_GMM_estimParams1.mod \
+	estimation/method_of_moments/AnScho/AnScho_GMM_estimParams2.mod \
+	estimation/method_of_moments/AnScho/AnScho_SMM_estimParams0.mod \
+	estimation/method_of_moments/AnScho/AnScho_SMM_estimParams1.mod \
+	estimation/method_of_moments/AnScho/AnScho_SMM_estimParams2.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 \
@@ -1275,6 +1281,8 @@ EXTRA_DIST = \
 	lmmcp/sw-common-header.inc \
 	lmmcp/sw-common-footer.inc \
 	estimation/tune_mh_jscale/fs2000.inc \
+	estimation/method_of_moments/AnScho/AnScho_MoM_common.inc \
+	estimation/method_of_moments/AnScho/AnScho_MoM_data_2.mat \
 	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 \

tests/estimation/method_of_moments/AnScho/AnScho_GMM_estimParams0.mod

+% Test estimated_params block using initial values without bounds
+
+@#define estimParams = 0
+@#define MoM_Method = "GMM"
+@#include "AnScho_MoM_common.inc"
\ No newline at end of file

tests/estimation/method_of_moments/AnScho/AnScho_GMM_estimParams1.mod

+% Test estimated_params block using initial values and bounds
+
+@#define estimParams = 1
+@#define MoM_Method = "GMM"
+@#include "AnScho_MoM_common.inc"
\ No newline at end of file

tests/estimation/method_of_moments/AnScho/AnScho_GMM_estimParams2.mod

+% Test estimated_params block using prior information
+
+@#define estimParams = 2
+@#define MoM_Method = "GMM"
+@#include "AnScho_MoM_common.inc"
\ No newline at end of file

tests/estimation/method_of_moments/AnScho/AnScho_MoM.mod → tests/estimation/method_of_moments/AnScho/AnScho_MoM_common.inc

@@ -20,13 +20,6 @@
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 % =========================================================================
 
-% Define testscenario
-@#define orderApp = 2
-@#define estimParams = 1
-
-% Note that we set the numerical optimization tolerance levels very large to speed up the testsuite
-@#define optimizer = 13
-
 var c p R g y z INFL INT YGR;
 varexo e_r e_g e_z;
 parameters tau nu kap cyst psi1 psi2 rhor rhog rhoz rrst pist gamst;
@@ -67,7 +60,7 @@ INT = pist+rrst+4*gamst+400*R;
 end;
 
 steady_state_model;
-  z = 0; p = 0; g = 0; r = 0; c = 0; y = 0;
+  z = 0; p = 0; g = 0; R = 0; c = 0; y = 0;
   YGR = gamst; INFL = pist; INT = pist + rrst + 4*gamst;
 end;
 
@@ -139,95 +132,93 @@ end;
 @#endif
 
 
-% Simulate data
-stoch_simul(order=@{orderApp},pruning,nodisplay,nomoments,periods=750,drop=500);
-save('AnScho_MoM_data_@{orderApp}.mat', options_.varobs{:} );
-pause(1);
-
+% % Simulate data
+% stoch_simul(order=2,pruning,nodisplay,nomoments,periods=750,drop=500);
+% save('AnScho_MoM_data_2.mat', options_.varobs{:} );
+% pause(1);
 
 %--------------------------------------------------------------------------
 % Method of Moments Estimation
 %--------------------------------------------------------------------------
 matched_moments;
+
+%first-order product moments
 YGR;
-INFL;
-INT;
+YGR(-1);   %redundant
+YGR(0)^1;  %redundant
+YGR^1;     %redundant
+YGR^1*INFL^0*INT^0; %redundant
+INFL(+2)^1;
+INT(-2);
+INT(2);    %redundant
+
 %second-order contemporenous product moments
-YGR*YGR;
-YGR*INFL;
-YGR*INT;
-INFL*INFL;
-INFL*INT;
-INT*INT;
+YGR^2;
+YGR(-1)^2; %redundant
+YGR*YGR;   %redundant
+YGR(1)*YGR(1);  %redundant
+INFL*YGR;
+YGR*INFL;  %redundant
+YGR(2)*INT(2);
+INT(2)*YGR(2); %redundant
+INT(2)^1*YGR(2)^1; %redundant
+YGR(2)^1*INT(2)^1; %redundant
+INFL(-2)^1*INFL(-2)^1;
+INFL(-1)*INT(-1);
+INT(0)^1*INT^1;
+INT(0)^1*INT^1*YGR(2)^0*INT(-2)^0*INFL^0; %redundant
+
 %second-order temporal product moments
 YGR*YGR(-1);
-INT*INT(-1);
-INFL*INFL(-1);
-end;
-
-% get indices in declaration order
-iYGR  = strmatch('YGR',  M_.endo_names,'exact');
-iINFL = strmatch('INFL', M_.endo_names,'exact');
-iINT  = strmatch('INT',  M_.endo_names,'exact');
-% first entry: number of variable in declaration order
-% second entry: lag
-% third entry: power
+YGR(3)*YGR(2); %redundant
+YGR(2)*YGR(3); %redundant
+YGR(-2)*YGR(-1); %redundant
+INT(3)^1*INT(-2)^1;
+YGR(5)^0*INFL*INFL(2)^1;
+YGR(5)^1*INFL(-3)^1;
+INFL(-3)^1*YGR(5)^1; %redundant
+INFL(3)*INT(2);
 
-matched_moments_ = {
-    %first-order product moments
-    [iYGR       ]  [0   ],  [1  ];
-    [iINFL      ]  [0   ],  [1  ];
-    [iINT       ]  [0   ],  [1  ];
-    %second-order contemporenous product moments
-    [iYGR  iYGR ]  [0  0],  [1 1];
-    [iYGR  iINFL]  [0  0],  [1 1];
-    [iYGR  iINT ]  [0  0],  [1 1];
-    [iINFL iINFL]  [0  0],  [1 1];
-    [iINFL iINT ]  [0  0],  [1 1];
-    [iINT  iINT ]  [0  0],  [1 1];
-    %second-order temporal product moments
-    [iYGR  iYGR ]  [0 -1],  [1 1];
-    %[iINT  iYGR ]  [0 -1],  [1 1];
-    %[iINFL iYGR ]  [0 -1],  [1 1];
-    %[iYGR  iINT ]  [0 -1],  [1 1];
-    [iINT  iINT ]  [0 -1],  [1 1];
-    %[iINFL iINT ]  [0 -1],  [1 1];
-    %[iYGR  iINFL]  [0 -1],  [1 1];
-    %[iINT  iINFL]  [0 -1],  [1 1];
-    [iINFL iINFL]  [0 -1],  [1 1];
-};
+@#ifdef Extended_Matched_Moments_Checks
+YGR^3;
+YGR(-3)^2*YGR(-3); %redundant
+YGR(-3)*YGR(-3)^2; %redundant
+YGR(0)^1*YGR(0)*YGR; %redundant
+INT(-2)^2*INFL(-1)^4;
+INFL(1)^4*INT(0)^2; %redundant
+INT(0)^2*INFL(1)^4; %redundant
+YGR^2*INFL(-3)^4*INT(5)^6;
+YGR^2*INT(5)^6*INFL(-3)^4;%redundant
+INT(5)^6*YGR^2*INFL(-3)^4;%redundant
+@#endif
+end;
 
-if ~isequal(M_.matched_moments,matched_moments_)
-    error('Translation to matched_moments-block failed')
-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
+      mom_method = @{MoM_Method}           % method of moments method; possible values: GMM|SMM
+    , datafile   = 'AnScho_MoM_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
+    , order = 2                           % 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
+    , verbose                           % display and store intermediate estimation results
     , 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
+        , 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
-        @#if mommethod == "GMM"
+@#if MoM_Method == "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
@@ -249,8 +240,12 @@ end
 
 % 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
+    @#ifdef NoEstim
+    , mode_compute = 0
+    @#else
+    , mode_compute = 13                   % specifies the optimizer for minimization of moments distance
     , additional_optimizer_steps = [1]    % vector of additional mode-finders run after mode_compute
+    @#endif
     % 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
@@ -259,7 +254,7 @@ end
     %           ,'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
+    %, 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
@@ -277,10 +272,8 @@ end
     % , 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                             % 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
-

tests/estimation/method_of_moments/AnScho/AnScho_SMM_estimParams0.mod

+% Test estimated_params block using initial values without bounds
+
+@#define estimParams = 0
+@#define MoM_Method = "SMM"
+@#include "AnScho_MoM_common.inc"
\ No newline at end of file

tests/estimation/method_of_moments/AnScho/AnScho_SMM_estimParams1.mod

+% Test estimated_params block using initial values and bounds
+
+@#define estimParams = 1
+@#define MoM_Method = "SMM"
+@#include "AnScho_MoM_common.inc"
\ No newline at end of file

tests/estimation/method_of_moments/AnScho/AnScho_SMM_estimParams2.mod

+% Test estimated_params block using prior information
+
+@#define estimParams = 2
+@#define MoM_Method = "SMM"
+@#include "AnScho_MoM_common.inc"
\ No newline at end of file

tests/estimation/method_of_moments/AnScho/AnScho_matched_moments.mod

+% Test translation from matched_moments block to M_.matched_moments
+
+@#define estimParams = 0
+@#define MoM_Method = "SMM"
+@#define NoEstim = 1
+@#define Extended_Matched_Moments_Checks
+@#include "AnScho_MoM_common.inc"
+
+% get indices in declaration order
+iYGR  = strmatch('YGR',  M_.endo_names,'exact');
+iINFL = strmatch('INFL', M_.endo_names,'exact');
+iINT  = strmatch('INT',  M_.endo_names,'exact');
+
+% M_.matched_moments has the following structure:
+%   - first entry: index number of variable in declaration order
+%   - second entry: lead or lag
+%   - third entry: power
+
+matched_moments_orig = {
+    
+    %first-order product moments
+    [iYGR       ]  [0 ],  [1];
+    [iYGR       ]  [-1],  [1];
+    [iYGR       ]  [0 ],  [1];
+    [iYGR       ]  [0 ],  [1];
+    [iYGR       ]  [0 ],  [1];
+    [iINFL      ]  [2 ],  [1];
+    [iINT       ]  [-2],  [1];
+    [iINT       ]  [2 ],  [1];
+    
+    %second-order contemporenous product moments
+    [iYGR       ]  [0     ],  [2  ];
+    [iYGR       ]  [-1    ],  [2  ];
+    [iYGR  iYGR ]  [0    0],  [1 1];
+    [iYGR  iYGR ]  [1    1],  [1 1];
+    [iYGR  iINFL]  [0    0],  [1 1];
+    [iYGR  iINFL]  [0    0],  [1 1];
+    [iINT  iYGR ]  [2    2],  [1 1];
+    [iINT  iYGR ]  [2    2],  [1 1];
+    [iINT  iYGR ]  [2    2],  [1 1];
+    [iINT  iYGR ]  [2    2],  [1 1];
+    [iINFL iINFL]  [-2  -2],  [1 1];
+    [iINFL iINT ]  [-1  -1],  [1 1];
+    [iINT  iINT ]  [0    0],  [1 1];
+    [iINT  iINT ]  [0    0],  [1 1];
+
+    %second-order temporal product moments
+    [iYGR  iYGR ]  [0  -1],  [1 1];
+    [iYGR  iYGR ]  [2   3],  [1 1];
+    [iYGR  iYGR ]  [2   3],  [1 1];
+    [iYGR  iYGR ]  [-1 -2],  [1 1];
+    [iINT  iINT ]  [-2  3],  [1 1];
+    [iINFL iINFL]  [0   2],  [1 1];
+    [iYGR  iINFL]  [5  -3],  [1 1];
+    [iYGR  iINFL]  [5  -3],  [1 1];
+    [iINT  iINFL]  [2   3],  [1 1];
+    
+    [iYGR            ]  [0      ],  [3     ];
+    [iYGR  iYGR      ]  [-3  -3 ],  [1  2  ];
+    [iYGR  iYGR      ]  [-3  -3 ],  [1  2  ];
+    [iYGR  iYGR  iYGR]  [0  0  0],  [1  1  1];
+    [iINT  iINFL     ]  [-2  -1 ],  [2  4  ];
+    [iINFL iINT      ]  [1   0 ],   [4  2  ];
+    [iINFL iINT      ]  [1   0 ],   [4  2  ];
+    [iYGR  iINFL iINT]  [0  -3  5],   [2  4  6];
+    [iINFL iYGR  iINT]  [-3  0  5],   [4  2  6];
+    [iINFL iYGR  iINT]  [-3  0  5],   [4  2  6];    
+};
+
+% Removed duplicate moment conditions
+matched_moments_no_duplicate= {
+    
+    %first-order product moments
+    [iYGR       ]  [0 ],  [1];
+%    [iYGR       ]  [-1],  [1];
+%    [iYGR       ]  [0 ],  [1];
+%    [iYGR       ]  [0 ],  [1];
+%    [iYGR       ]  [0 ],  [1];
+    [iINFL      ]  [2 ],  [1];
+    [iINT       ]  [-2],  [1];
+%    [iINT       ]  [2 ],  [1];
+    
+    %second-order contemporenous product moments
+    [iYGR       ]  [0     ],  [2  ];
+%    [iYGR       ]  [-1    ],  [2  ];
+%    [iYGR  iYGR ]  [0    0],  [1 1];
+%    [iYGR  iYGR ]  [1    1],  [1 1];
+    [iYGR  iINFL]  [0    0],  [1 1];
+%    [iYGR  iINFL]  [0    0],  [1 1];
+    [iINT  iYGR ]  [2    2],  [1 1];
+%    [iINT  iYGR ]  [2    2],  [1 1];
+%    [iINT  iYGR ]  [2    2],  [1 1];
+%    [iINT  iYGR ]  [2    2],  [1 1];
+    [iINFL iINFL]  [-2  -2],  [1 1];
+    [iINFL iINT ]  [-1  -1],  [1 1];
+    [iINT  iINT ]  [0    0],  [1 1];
+%    [iINT  iINT ]  [0    0],  [1 1];
+
+    %second-order temporal product moments
+    [iYGR  iYGR ]  [0  -1],  [1 1];
+%    [iYGR  iYGR ]  [2   3],  [1 1];
+%    [iYGR  iYGR ]  [2   3],  [1 1];
+%    [iYGR  iYGR ]  [-1 -2],  [1 1];
+    [iINT  iINT ]  [-2  3],  [1 1];
+    [iINFL iINFL]  [0   2],  [1 1];
+    [iYGR  iINFL]  [5  -3],  [1 1];
+%    [iYGR  iINFL]  [5  -3],  [1 1];
+    [iINT  iINFL]  [2   3],  [1 1];
+    
+    [iYGR            ]  [0      ],  [3     ];
+%    [iYGR  iYGR      ]  [-3  -3 ],  [1  2  ];
+%    [iYGR  iYGR      ]  [-3  -3 ],  [1  2  ];
+%    [iYGR  iYGR  iYGR]  [0  0  0],  [1  1  1];
+    [iINT  iINFL     ]  [-2  -1 ],  [2  4  ];
+%    [iINFL iINT      ]  [1   0 ],   [4  2  ];
+%    [iINFL iINT      ]  [1   0 ],   [4  2  ];
+    [iYGR  iINFL iINT]  [0  -3  5],   [2  4  6];
+%    [iINFL iYGR  iINT]  [-3  0  5],   [4  2  6];
+%    [iINFL iYGR  iINT]  [-3  0  5],   [4  2  6];    
+};
+
+
+
+if ~isequal(M_.matched_moments_orig,matched_moments_orig)
+    error('Translation to matched_moments-block failed!')
+else
+    fprintf('Translation to matched_moments-block successful!\n\n')
+end
+
+if ~isequal(M_.matched_moments,matched_moments_no_duplicate)
+    error('Removal of duplicate moment conditions failed!')
+else
+    fprintf('Removal of duplicate moment conditions was successful!\n\n')
+end
\ No newline at end of file