From 562a9c737fd10a91679ff86a55c3da2cc8c0969c Mon Sep 17 00:00:00 2001
From: Willi Mutschler <willi@mutschler.eu>
Date: Wed, 6 Jan 2021 14:03:51 +0100
Subject: [PATCH] MoM: Improve testsuite

- add Andreasen, Fernandez-Villaverde, Rubio-Ramirez (2017) test models
- move models to dedicated folders
- add `make m/method_of_moments` and `make o/method_of_moments` commands to run testsuite only for method of moments
---
 tests/.gitignore                              |   6 +-
 tests/Makefile.am                             |  21 +-
 .../method_of_moments/AFVRR/AFVRR_M0.mod      | 299 ++++++++++
 .../method_of_moments/AFVRR/AFVRR_MFB.mod     | 300 ++++++++++
 .../method_of_moments/AFVRR/AFVRR_MFB_RRA.mod | 299 ++++++++++
 .../method_of_moments/AFVRR/AFVRR_common.inc  | 540 ++++++++++++++++++
 .../method_of_moments/AFVRR/AFVRR_data.mat    | Bin 0 -> 10941 bytes
 .../AFVRR/AFVRR_steady_helper.m               |  80 +++
 .../{ => AnScho}/AnScho_MoM.mod               |   0
 .../{ => RBC}/RBC_Andreasen_Data_2.mat        | Bin
 .../{ => RBC}/RBC_MoM_Andreasen.mod           |   0
 .../{ => RBC}/RBC_MoM_SMM_ME.mod              |   0
 .../{ => RBC}/RBC_MoM_common.inc              |   0
 .../{ => RBC}/RBC_MoM_prefilter.mod           |   0
 .../{ => RBC}/RBC_MoM_steady_helper.m         |   0
 .../method_of_moments/RBC_MoM_steadystate.m   |  74 ---
 16 files changed, 1537 insertions(+), 82 deletions(-)
 create mode 100644 tests/estimation/method_of_moments/AFVRR/AFVRR_M0.mod
 create mode 100644 tests/estimation/method_of_moments/AFVRR/AFVRR_MFB.mod
 create mode 100644 tests/estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod
 create mode 100644 tests/estimation/method_of_moments/AFVRR/AFVRR_common.inc
 create mode 100644 tests/estimation/method_of_moments/AFVRR/AFVRR_data.mat
 create mode 100644 tests/estimation/method_of_moments/AFVRR/AFVRR_steady_helper.m
 rename tests/estimation/method_of_moments/{ => AnScho}/AnScho_MoM.mod (100%)
 rename tests/estimation/method_of_moments/{ => RBC}/RBC_Andreasen_Data_2.mat (100%)
 rename tests/estimation/method_of_moments/{ => RBC}/RBC_MoM_Andreasen.mod (100%)
 rename tests/estimation/method_of_moments/{ => RBC}/RBC_MoM_SMM_ME.mod (100%)
 rename tests/estimation/method_of_moments/{ => RBC}/RBC_MoM_common.inc (100%)
 rename tests/estimation/method_of_moments/{ => RBC}/RBC_MoM_prefilter.mod (100%)
 rename tests/estimation/method_of_moments/{ => RBC}/RBC_MoM_steady_helper.m (100%)
 delete mode 100644 tests/estimation/method_of_moments/RBC_MoM_steadystate.m

diff --git a/tests/.gitignore b/tests/.gitignore
index 722a27d5d7..275d8736e9 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 fd050916a1..1b7929a582 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -50,10 +50,13 @@ 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/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 +838,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 +991,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 0000000000..8e51ac5136
--- /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 <http://www.gnu.org/licenses/>.
+% =========================================================================
+
+% 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 0000000000..450739ad3b
--- /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 <http://www.gnu.org/licenses/>.
+% =========================================================================
+
+% 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 0000000000..9c069d3a3d
--- /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 <http://www.gnu.org/licenses/>.
+% =========================================================================
+
+% 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 0000000000..76aea9e0b0
--- /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 <http://www.gnu.org/licenses/>.
+% =========================================================================
+
+%--------------------------------------------------------------------------
+% 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 0000000000000000000000000000000000000000..f606b2109e2bf61024b860d13000842c27f3d290
GIT binary patch
literal 10941
zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2cQgHY2i*PhE(NS<NN=+<DO;O0t
zvr-68O;PYl%u`S>Q!q5JGB&X?HCHe)Ffvpi5-`93qo*%F0|UcORtAQOIl2?MO|F$o
z9FIS_JHjn>qn=aAl^AU`-_DZ(>kQX+FG`s<C1i@rdX1A_Cyr!iP4(2At2lk%$~|ki
zgJ<bI>JWJu^5}2uzo~!Auf3o1yr%xV<@xtNZT7$a-0GqvsXWzV$^Vs4>i@o9QG06j
zk$|5f^E)4k?lC=KTC~$8cF%uZv9;S{7EIkcSHb?n^#{KbS@zAgoyT>NWu@?1FRtgp
ze+&zb+_^X<q3p&yh1&^C?cGa_&hs=){IE&dt(!CF{*~#k+k`y+7Fk7eAF+Bk>)G3+
zCXK13m%6H&o*mj)W@XFreOIcuT)q3mfG5E(pLz+-KXh&X6PE4zb!8P_xtT7R-C3%)
z@PxF~nh+C#8UAmIUq?AlNZsoHZ7ci1V{UizPS5u8G0V8PLU8MgRg>@TJnbd1UUOeS
zbyR27<v-=y&uJM>2-@)Qq^)Y2J4>wny;Vn~r)Nqozda?w{&aD8x16K%-jLnLIv%WD
zeeB$^k4pQC|D-BwHx)1YEMsKK+4l6|y?tdYimAnQMP`Z~*;9L~{(5oxx_*$oY1PC#
zwO#72GH2HDc~Z}(IwV|~tYTk2;XvlOCz)o2EZ_6*e$CTV^>`PO|1ZHq@#E^^+eYn&
zZag{|dRJ6L?QJP{@@CCNBG&sh?{U1iWc~YVSA8}bcQ4CYzobBZ>q&**(Oh>I1!}Hd
zA`*S|wvY2|&Utz}e=jlVZ3wFhHwzJ#h}&(VkgKL)zUuPKD~=Z)Zd<38Cw6qE&BI@f
z9$f-&e^1Iz*T^a1U6$27=ia2=_aeR`t3=n;SviYbHd6lQwLQ_pR>A)3fq98>GY)qB
zdoL+8souoEKiOj9{kQUyjwd|HG}GTO>7n>e=dUly*fyozT)gru&*SgsWpgqTUu<0_
z9xrATGNFTSYeC!8?z__Q{Q*9aA1<sa);s)L{?(TG>r-~<Y6nhHHFD8gz4v?E8T}Tk
z7gOd7zt)}hZQl0e7i+G@m^iuam%Art{f_tObBT;qJDd%Zo|-<3vpE{w{^8{0O&W1i
z|44p#rskS3y{7a+^2FqTsJHtjCuUr!`06(EgtW(|zFS8U#2()W7SsJ?Vxy)JBe5oQ
zRrpVxlpFK^1eDGZbDnbf!ujd0=N^jdnfWb}*LRxZtCwYP<67G3vr!tq9&OG2U96Gz
ze2RM8^&lJL>(VFJ3V(n0|HLC{R@cOff*H08Z`3KM<<52Zcu`mH@hhITJ<9(+a{0Gy
znryGSwl~$H^5pjy1saEE{JXbc=1m@1$=qPBZQY+O6EA+zReB$mn^z~$ym$7)xS1B|
z_LI+ci{58XIo7`}=Foz|Pp9}Q^<5jQ&M2}^e7)!S_p@5IB9pgFH!crlJAdWaUB+KL
zV$9xWF7ZmR+1*+b`$Fg6<>guNJ?$^L*B`rbhlOoEd%M-=)e_fssb6F7_-$<UtL%dp
zyL10PwQdXEz6`a8b=sCo82??oy{#Z6V)oIemW3~-oVoJS`{9di(j}%o4<cfI>^ZMD
zYq9BOm#^nE``%yv@<mC$IK<*jUHZjyiT0O&hCHh7J8rkJs=tKwPto>!7piwP804RQ
zGxr-)&slbzbvmL8uBXczsrxycJ<PA3{Eb^9y=(hFF2;#h``TVbGIh-sizz%O{2}$F
z_3g*$tSWQLH_p+{Dwy|i+KFVQ89x_X@7&9@XVrgCyB)p`=XjbQED~Y-8GP=`N)e75
zA!+hE%)=eRdruYo2xOEGO}lDp#qD9Av(d1o$-@8h)1wo!n5G;)e`)nCi<HgARj>Vv
zYeW}+EpnGSzH!NuJ9qBN2(Am97PU`i*`;45yIE!afBD`vbw<&b6>c_4i82eX8PEOm
z$ik=VX|>u~@!Dv)70cw0Z1m31zZ7+_Fd%+rTSVV|ndqAvYh|8I&wDjxS(lx8?&Wtv
zzs+qocz-|f<i=jZ{8fD`0z2a)!&M9a?6C~jx1M)6HzVoe<+n25s~=B0!~bx@C*v^f
zZ;w4D_J1rccyi(MJG-JQEstxP_dah?YpZ8mIXiPVr~rQjDZnRJnoQwlINI)L>3TjY
zYKoVw=0~T=RqM=rR)}8;xGAI>-63*vT4DfO@@zx5v^UXP0&05II>Q(By$x(R6x`Xh
zsLE=6^||+Vs?YEI9~UXqX%gS%GU?aLCI9~}SE>Ix%koIw$IbEu---=Bmfp~{oj(6u
zZ$_%<t;K$+0kx|=R!>cPeDcwbjR7lVI!t3WNQfjSy?p)q^^}R-t-N!V9i6n1qvqbn
zcQ20e9BDC`9R2-?)7n>8=W6G53d~<MEwXyEjQ-rXsMYiJ8fO_C>|Yq!wKOVa{+p+E
zE*s`tv?<C^b?eukTa$RwdBs()-@Dzr1W)o_*UH=7yWL%^^2_GWhWpM=dM=*c@_g%}
zYjNpv$C7U8yR7->;=SkWInOMero5<j?^4c!PcxQ(oA-9&g9|sB*KCw9*(Q^xle@Q3
z>SJZ@)VZw@9+$SM>^Xh$v1(4uBGH~{uHh*=r&<0Ge0u5j(@7ka_YbQ}9XWLIjj={b
zjn|>DUvHZG7Fwk6eAhbS%dWh=BTiz^>{ih~y|-?bzHFIqw=iqdCYjJAvm)b_FMr(7
z_MiL7-|fyf<wLpNT>+WZi`s6UZ@+up|JW2((H-LFzur%a<LJEbyE?o->>$(jKgnND
zJ$x*xpYyp(>DH^T&hlfw4*L2UE%|HJy~1|$jnJ!34k-a1`$Xd!KLpwd@0}&plRBGE
zzQ!iXaYfAQ<UM6R?bGf`ul~14=*!)iX7~SlHH3W@_hP#nvu4lli&rWZF$r}=f62LG
znKA49#kEgT4u?vWPn9nIxFz(>qW>#p_GJ4N$j6Cx_x`P7bd$=F_^>}~huhP|cKbW^
zI8q*%+&`w@9n|@`I^8AOkL^%_%IQ4?{zmJ5to?r>RHX5oy>pM9Zm;62FTT<3-!?gh
zzX)7<+gh;x_RF$hk;xqu#>+ELZIf!+c|hUeFGIbOna9GrQh!9;d_3u=UpM2ViJsAu
zOC<vTzu0mw{dcj!uGuH{=N?zT_Fnn4)laoU3odOqkb2OvM(@Wrx5?9(CLVU`+`ck(
z!#;7#2|34DXWvr(pL(^l=l`@TS`+k}kExgQm$xcaD4lxnubbzM{ELuNVY&y-n!m~n
zl#VF)tyo@qIMS{zUWI4XS1s53J|6$x*zu`df3wb3Plx^KqRDqlw%UGlDa$#iJHKQ_
zR$kPW#Z%hTzi-)}b4uIwW0}{AH?J%-y3QTjp#F$s-HQG*cjvT6pAS)6|Mla89}hjf
z5_U-@t_+x_aa4~nSn<G-Z;o7Rrg6+VxjN#-sag|J(J9T|Cr(8w*UIImtJ|#nD7av<
zi0l2wJ2q^neY3TtiZduN^mNbfh>BbOT-Pd}bUe(FoEflB_tF%exy%24?wDlEQ}o_J
zDo}Vzo<;C6mX+<xEuzaLgi?HYC%w>``2N1&lFQLOdB0A`&sNrJzEUt(aj!P_Blh6q
zcV#Oxe0~+*7r83Ee$v&$GG~k~t_zA4snz9{UZel=%Dkk;n%z~|&vVQ^*_Oy^S<NW;
zTjb*Yu~O*LwkxMkpE>XHV1L=V#rjPgK~GXUZ#C<Pa81cDmkK_Zzw5q0wFaX!``Km2
z&&!f8{Ml$at9|BE!Rl9A7jF!x`v1J)gGz*s<f{!;#~M?O@-tfhN9ypeT>H;`odmz>
z7s>wGO%F>SePx`kz?Ua2)y{KQqtVPbR^Zko?qyw}AFp@Stg%@2v@EWbS+$b0`1wI6
ziQ<>-a_deszVP+Acl<2Vi^%PDnT47Sbqr_sJ+zR%%rCO3&S!I?y=mpe*FjQWjBjUV
z`8jgLXYV}ytWd#+S@P*JU&VlZYZizZeh+@$5gX!ZA0JZsw4VEh&+$KpHhOJ7;WFbo
z|EXh__b8^!oPR}sf%DSsMbE2~`OaNXXIuHM(!lff{727xkH7l<$gC_iQX;d@{_DZ1
zjq5IEY|Oe_vPFB9ch#%s3ldap-a8g>?TYMU6p7WFcHOt+^$gWU6|2+MzaB=<J>Ict
z+uL%vmY$CHrAIe@h<~=q=ubG)joF7v*0&ZYe3-~DU3R|Zh}FR#e5W{juC0;H*u>3#
zw0WJR+XaI-r%=684~2is*F#Uw1C{uHAtnBlS_@$gnd8fy;;Jl*id*N1dTYhZS-JbZ
zURj*i(&{%)mhO@{$t$={INMWYY1`A(E0ZRg-0qQ_T_=%u_0+qN;PTi<R{G)E&u7HN
zgrC3rz0&Ud_usYi-rv99%EFYmPp~oJ%RSEjKdUAG-M@V=@q_*Q-`n=Gm;Y9M&u36M
z-SoZNmE5&)iPpCs=CrA||Jy7v=b*|QpGnSwYajev7iMOYccx}jLvZWDRhMp_5|XQV
zV*KT=q$9t8x`^Z9_N5;lPG0<%zgOhLi%o0e5An}6s5O1}Yjda28=FI4I&1>i^QUFX
zswDpu*na7j(e>@@Y{|QxzI_zN5hc~<Cx2QY`^cTSTb!BXcCt2aUgNi6OU>8N7m12S
zW!c)l9{DA8_W7Ojcz=ZBx^tunyQSFXR(<Iex&j%O4t|te+aRv}UdFgL&|zQUrr3#7
zMciIjEa^Jq^~ctI_xWenQavp8diY#0clPpln>u|_na9F*d5NX6%t~|f-j}@6N(s2f
z6L;rI-~FmKxrn`sc*-TON-XQ$@~H2(w~tvhhr_CwH<EX&vghB8Dfn|cRA|A37PlgQ
z)qq9L_8MCR9j^-*U!AMftQOAR#<-Q;=hunf`jy;@YF!^wiXXWa{IM#V`Kge<=G6b>
zvX?zVN)a#B_Qu6LIDT%))1{)1l{dTG%`TD^k}sOIEIsPuhm)1Lc}JDbNcWk2S&@Gz
zXf@ZxZO*kGk{cd~_hoN*@?Jvz_0y!1OY1*H+Fjs@ZPhDtV-Za7xbpLXm8wVXjh)uf
z(hkLMO468`*>}h6VAWosb)v|uCh_AQW;^M9Q$F_8GK;CEx<zjdan@XIXn1ff<At!;
z(mA0_n(kW-Cmysvxbvk$Y*2#q3o%P+1D4(kS8E$hW}IrQEoT#*c=GBFwfRAx&EhyV
z_GSjR>@7WE@J(*9TIa;AM#pDctUX$|`@y%lT6#yGTKY-9f7$d(D>K6{!0^NG=PkK%
z35R#E2iDfdn$2>T-Rr(+`r_HIzo$u+vj2OuceCB;@Xlj*V}GuSaCyHZ{?PGHH+I~*
zuNAz0a>Fs@Hy#(wV?Xrd7>HlKS}=Fk*{k30vpn?8wg2v8{pDWXguA-oO$}CuOutV*
zaJg~9$*63_kkHeo=45$0yv-@Rw0<`K(w9^B&F21Lp%;6-+&q-U_$a5e{fmxc3vXy%
zTQQ~gd&<r;E1PZBd==g7FEY1vR^|NP>hq>1t>dliNYib)AJo-mk-}Ls^VOt%3v}4e
z%cgEWqdf1BKmV1&Cb<sB&|fa=q&%3J+V@=w3Fj=fkKY)ce5{#c)_wb)HxfA~e0QAN
zlGSjg^4;d^$5c<0UJ;*B9(MScJhRc4y}2uP7Wv8V4{O+4%jTTTBKzaA)xBojNH3Ls
znaw`C;yc!vv$%frisvb_+~JdU!i8D1GjP$G7GIs{eC2k%ip9_8Zp>WYy0#{~@(Q!%
z1^1A5-hJ~Mj;?)u^i?6}8taGGC1)!z81C$OI{h~9E6-Uai>2j_f;#F7InO?by_R1y
z`y}gj?Rl-QJS1&m50%8r-eu|B<@K(l*7BC-4tKZBs~^^M?E1aXiRbGT^PZLA@BbzI
zd3V6_arT8@or!$sS93bGJ<`^jaiw{0Yu!v`M;^cD?kR_t30YJuyxp$&zD4WQ{B5iU
zWio=R_S(hhu?J_T%g^VPH4r-xs9xQrP@A{v$95yj6T6baKlr#N*uRO1@y=rVyYjB|
zw>P;>*6W+>PR7pR71Xz$nY#Vfg8R?EUbYktno;kOxIg)c^^LfqsXNkE+Dq6c$?g37
z)b&r^m!N5~n*~dSL*h1n<2vWBCy{eOqhQlgz4->w&HJr>#G1_4JgWUpKQuct=S%F1
z{S^!Q`Zpcd-}iwxS@D1NwcPK=#5(SX7hEsgCfgSJ{l#-nrwL{?dlZcy=-hf<yDjcZ
z5TkedEc-dzCLEAozE#flHEY6t`}@}OB6$rj#zmi9)^gy^t|zaejqO<6RNY$Qz1s?A
zExvqmquPwk{|}zC|Go3ab%)9a#ft>w{najRyC;w+t@^wzIYRW_wT$WCe-$MLdEB%M
zF5@yc)G+lqUDXrwdt1`!wQhgJuXlYE1XThrpp`&%h%gKD(c|9cl^j!7yb$j)@D`Z5
zVM@d7Cp<3Q?yVh6dJ3Nym73J~A5NH_Y?GAC&-Ugdo3w?*85T~VfFnoTnN9z4TyROy
z;9;G<MA@eF{Li`7?=8Q-&vsT2=-}pPVa=cOVg3GzN9+F{vQhj0aap~iZy%ROhqceF
z^V^@fyElcL`XkQSaqsDK(Y;-rbEMWko#-!V@oJ)xiTOm{u7EwupS+U%Q4w$dTks{1
z2QxR@naM$k&(?kap0c^%!b$C~J2d1w3w_m@mAp3WTl;63tB;zK(CkxZKS_3^ZIquk
zZ_(pl>yNSb=YCE|u5x<3X+rN3b*aAX(TsMN_6rLIdbfw33CJj1BETb8zh<dFOVJy1
zmqk$>i?<Z+UFF76`QV)ABs=l8xr;T6>o*<P9+maU#J@@TX8e0CP2t1pW@~3BX*TE0
z++F#9s#nMd&UWs*sVD5j_bq(OZrRt@BKzmcJ|WjiyM0$2RO@6uPFV45)y8dq{{4z@
z6pcC4*_>J$x}!hvqvtA_)rYEX*K|yMq308(V`H-a|Gv=oU-pSC>+0CQhIy5P=Cnz&
zJwacV+-Qren)XP0{+i(Vmr^2({=G9WZtgDKW&M`#aZ<)Rw^w{6oeQj{-%DP*&aL+S
z-498T%2(!|s(yWg*J5e3p5oVx1j)h`$14M+)8w;Hb}gF!anFfQm$uvKdQ@yNb~%)H
zjD20Cm;G<IFZ~<*#A8>w?GM?mpR<i|t!LymA9m9(Yg`*$Wh}ey{@bzhx?fM%TK`D~
zYn<MSJ)ahKylRbkFZVm19EW=EMYTCoq-I)A+P!>LNX56KyS}g2x8`;hsr`1Nz&2|c
zYv}%OnH&DTJ(O~vJ!{I*t1aj9@>l}xOmb4Mcif%+cAs?G7Ki;Cv$I#!8nf)?z9~I<
zw!rrNYic)#OSi0#lk2)*qWG7ic)s`3N{M%}HpvG0ZXr<*g}**Z`!WB2*SC*8wkG_s
z^Yd4)IL>fXV4B^srmgocO}u@VYeMzb{950sPL>s)gu9oQ7AU7)O&2>=7@cbvepX=d
zl~SXpw<fS0x0+?TrA*Z-GVWLF!PFU^rLjeKT=<^3_S^o84azt!XRF%B`S4-h?*Cy?
znq~Db!?agSc=((<S@Zczrj#%B`%1JM`x`%N%1r9z*s@nT$I34x;nU2miKX5;eI-BM
zC(h(_h>Df|Hhm-0lfO@s^QM{Jxc%esu5Z(mQvUdbKV3Ir;mcE8tETB6t2W%_Jv-@0
z!OwHIUK+acMHS@jE9+(A+&7W=ck~jEQl8|g%fk<so--~<a^ZC4_&QyfQ)U9&_N!$z
z8cp3x%L;<$bH%+ip82)@bi%VQktZj#3SCYY-~Ri>)Qo!DRXcJ{vrnA4)o^PU`^5SA
zAIz*3^L9+$ossSGL8@et&K^~ly7ecXE<aE-XW8oJNqL_GUp=^%zVP|O^8MW3_f|*D
z;ObmnwWE0AtAu;ex99S&cFk7)_g{>i|Ag-O=DGDQXNrVRpYF|BlVfXqDO&Z=n%N<j
zjwC)zyElzf+MREk|9!n5kHy4y2`M+e2=3YRf2OH$%Fz#U;;~cz9&6pWB*EKeKEG7|
zH^r@!IbZK;v5{LUbAFAr^7UPt9`1hjX=}`C8JY9Fa}s=>DO8``%ch^oXy#$g_ExL&
z$coRsPl9Iq_}n@-jrE9eg5YHx^Pe;L+)Dl=g~`gkmVGzDmw5+AeZ{#4;)jo2Ey+s?
zkYX==Yo&bAnJF&y%etGQK6?ywe#~C<?9rPoKGIRg9l!gYmvM}=Ty*cF+V9Vi74H+T
zeo8gs;;Y?x`Cr3BhAmy}ou#$Q4Oa?&l8Y5eN)6!uDjdjJ@@dh^-S>11mLK9yVb^{z
zV{`3<O9K7k+Ny4B)BCU0y<Pc7Qmtv<_uLvgv;Wa0%K1f;8LU^79-f)EG{bNA41*vA
zrt3e>mYrv|*|JWp&(eJAqtFd6znsjptgF78m~H)E*m!-MwCr}##rw^lHfNhM?b)(8
z^3IO1#)%^P57gFK*WLW{T&8@FchkkD--`}wyI9=%dU>N8r}%fHL*cP!0uoLhk)0L3
zjpaYX`QNwPK;``}M$kBYYe=+)OQ6VJqhrSFg(kK96g(MS`S(P-fbXfA2|-?Vd*pu#
z>|1i8U7%C%>|VPkqLMu7z8;g@Pb#ZEWmo*UKw{>}HI`fR%`eXh%f4Ux?Oyi$+P?J)
z9RH#mH~eMgf5`vV-kkmZ23_&vpSG-hwCdTZS--9=OKYEhw=Ahx&sHxo!r0@C&hPDa
z$~}G3w~EYB_tJTrr+IjivSAH7`=Rq%udCxqg@n&Ed%FA1;g5Bh%hpm*Cq3oD^s3C1
zxSQ`={l9yine8FJ>$aKE>Ex}ZH9faC=5K%fv-sQD6E=L=ysueLStYpdJ+a^Nmub-V
z<dA|@!oMdQe=RJ0yhr}?qLt5*+>?@K#IHzeP~O*ketzuJ{BZ5xEPp%ve>=T)u=^rh
z8PhGz>?fxpCai1PAvJUB)^Z<rxzB$t*oI108|*n``R%ZLDZ@nh%9V4!Ew?UTwft$@
zYgQ%xU*hkhLMDDXX!(U9(EQch;@MZk_MJFad!c%5sXjN~l)07q`@?;o8eMO<YJBm&
zzL>|>f?Z**w$$`<g3p9jq%GIT`>B#v`MDzNJ9p`#^DpK{%{9>zs4P5xpm6`UBdcfL
zbFTY;Wx}fiwnuk1D9gDX=`=Vi@>H+*-p(EM)|b|=VPBEP?j-&wXl{wzftN9>_NqIW
z)_<A1ul`HTi{Gx#+g`cYzpPxjclmxP<+wsn3YgCVN&&7@xC25$C60%)SNe3V5MCtN
zp`_TsrKQ01pF?Ge2uFxgh}<kUMox{XZA_2E^ecH)-lsO+-MxFm>35%uCQM(Esg&5&
zp_XhX(HNqlbl}ar(@)FN-_6;Ze!hHPaozv_-(P;Ju*vn6^Pk+a=jBA3;v4Tjf3fba
zJy}2RVBl}%wG}q9zskk;p7?)Y8vDGJJ2|erig{FPE`Mk$Z}6XsH6qq)e(JvWRGhcq
z%|X_?Oe{L(b&{XF6?sg`UPozhU!1V&r_=qtEbC6j+s$3V5+SZP=c>hWk(!J1R$mim
zwb=1w;r&NejCuJrdE37p4%+jq|J&CioD=gzlQVNd+7w>LzCNp|`0M3a^C{B~zAJOT
zaV%I#WS&v3wX#s6lGzqRD{mg{+?z)KnS*%F@80ll&dwR##TmaWj$5piG1v7_6Z~{;
zV!Kmc(L2G&O(K(|x0ms+om})Wv)`>}m+}*X?dSKDU5f2fd|!3=$ib{Uk)CMVl~GIE
zZk6qM+u*{u{>q&Xt(;sclC6JU{&y6e^H@Ia-^b-#Q{L9p2L|jGD7XIq?dRf{uFJ01
zUxtW&5YFGVFLgP?RAI(M+5giH_2mdhuxS@(1x{ZrF*zkl^JvrQ=dXiJMLk|xpXu%h
zPS&4kw~l4<zlmH=eb;vT)@B{J{(h<P3J2$Jvy1`?3x682F4dM!;kf-bHZ*sG%rzCA
zg<Sp3EGAiodPWXMH12KxHS1m~XV6&=%dW?>4F5+<ubj-mw#=@$Wu{Qc=eqf8C#$gQ
z-Cb>2>%=%cwV6M6nwCMTz}3wKOQaTB$1VNovijhPe|z{oXKuZ+?dzMm(u|#2H==&8
z;%r$|{+G2T<#%)OgU4cF;ieLI&iWS2m5|H1d8f~L-F){w9yzR$SGO*beeyc}epuN2
zv=05s3)Wi{thz1Ma^C#o@_%-();062JIXgqvTjsgS@`PQzvOvsTRAR%zP)B&P<PJ_
z&)@5<UOR>qOJyYNyppVNB}wP|KCbtMZ2`(goheIp3eT*W%8{l1U`j`K(P_EMQ<^W>
z-i+q7(&(PA(fHWhB~kO9#C|#JTl&ov2WsC+Klbu#e_0(?RG0YSR(SKli%VGw&-!*t
znTYGhmfv3@cJ-{YD_`~%gDdJSp?54EH@~b3bZoyF(D^6!=%rKEMH1Rd{OLTW6P$Oy
zz0M=KukQ2_9kW)a`}Ru`T5cS$aoUsZ?lC<fPfYav-nA31F6Wr{=}q(fc9BJobr_>~
z-u#>*)A}_~u;cckDzC!Fo1^`mb+rsG?NyFA-+#!sx0gFdMr-Gv9^a@X1<|D?W%IU8
zS-pSe*O`aaH)PIlvA!=Yp2Wbl>ajA{dVQ6IC+`;@*wXI7b$Yk*x9rk?^5-8|c{B!e
z6lotmeev!5IU83W*tXTBvFJjYNbDzp!n*zYu0QUYwpc{Zk9qGVPRVd_<F9tba>t&<
zI-PmF-A7+ie8az!iM6x)L|yJR@$8wqY3{xvwm#X)b{*E{<R1QIGCn^F6k?5S<~<cX
z|80N8k8NEq#U>t@r#^G;`4-O1)wMenwx93TOnAMv<+XT|^Q!bf)2B1u9ari*?DO|s
zB9CWb>AMpWPCVSp0y#EXJ;=EopT_!wUpM*A15h#b5ZW=}^tsF~a*VskbB_bNwljB;
z#-b$G8Nm;o0}rlDU1V_Xp=rW{8M%zc?7`|jl1zaPHY-G4xboks|8L8%HXvbzkG7n!
zh+tyU)Ca9aE1a*VR-e?rUw!hlcig?{=Dkh}Z$t!Yg#FigYX8?<=l^<XpBaDqw{~ur
z`-bJ(Ce<k?m+d)xrNsKDW9Cs0(=@Z|drMEPcKh`2K<CaBQ;Y5L+_FuT5B?-AnK^sQ
zkHx#z{|qxdu}##+Plta#>xmVi>=hNNyYELVe5yE6%c4<Ej_X<J=k@9z3!fgGaeZ;M
zO~sqNmVXX5@>fMYE8qLB=EaT8I!Cg9X1vLH7rA3%zwZ8T-~alE9(WmfZenq#(iXLg
z{PurhTD}G^vI&n}<G9gy>NNL>53=qQL{6B^ach-&)c4IyXXg0rXy-V3!uHn7$TU62
zTi>HEt4uU{vq`2%a`oh-nL5XBw}*GPvc3$qx$YgfsXY1oyC6%Q*`>9SKZ|&>&c{BU
z8~&j5-&~P)L8eot>CxA3XbV~{@UB%-4){>9?wWPek+ad+?yu`ZnvQOIziOA>@vXk?
zZN+v|AOBjoGX9-~&x2#4_xAABTW*S~>W@slzQ{JM_o1OXpY_xu5^{>{Kc@%mf62{y
zPwSQV`x)yFbVchKPO9cwu`EXKcbQn@Klw{-TINCmk^L35!X626Id^4_Iji=VHTMU3
zdES_+7ku9!xbd&z0v$`I&Kdq3fo;;AD?&YkpUbIiEO<HT#T7YGH|}4uySY7cDyLLz
z;Z#Xm&#~)(uKK|}ucsF#`*d0zXJ?K7A+#c8>r|bLkVThvKYWpYU&}CQnwFBER^sD&
z>nnjSN;>C1taRPN<&q{j-6veRS4`c)H8^>De*Envi>2l~lC!9Lxq$8Ngs*)kuO_Uo
zEAE)R>Zs(b(D2YOTeoJv$a{hDEx8Boo{nF;iv6kLU-A51eA5>n`D_`juXKCs;)V5|
z!pHKa?TNo|y`_28x|XMAx(}k(?bN=rg1xxv=i@(_Re!`PKja(_v;9`O{pR;~*AlO-
z|K<3)ns)(j#Kr7qbKO^H*V<pNkgN5%{rp~-)swW#Z)b#;ZZLgyT4sNq%iHVIcerYQ
zd7!#8=En0o9SdbImmDp$`jc{8teVsJ^DTCnb!)R3CQD5)ntY$(_Jo!CejB$)r2hH)
zlRI2!-<f~BUu{*KRz5w=eKu6%6OUc&-kC0^BVxWq2zM%O7H{LZ{BPf;mhA2&r5(2W
zo|`UrZx**qTr+`%H)ig!+H6&+Gj@k0dZ%Y4tdQ8QH*2A5jz|AL=Lurl^;AOarmlYQ
zrTEO{Fmcm+GdBK-S>W_z>0Iql=PZ*p@2t;8S9kc-{5P_VYrYe{CFPjA=saV)W$Wca
z4sKpqdz$TXEzk37J9)z<L`Wvd&U?4?$O{#AUyXK_-m^Skqg@wS6fL@J-+f^6%iBll
zEwvuLx^VE^G_4IkPd!tMcb`08_l*&kWl+c)oh|X7Hi`YRlx-B1Wf9W(w8+@ku6V<{
zDt*xj3Y#zQZd3Qz^Rd$IVceo4>G^NtpIT^sD5=<4zc{e($OGneieW|;^^&jS=Pf+&
z=-`is`O}0rI$j!`TPAd&#qxLL1`&@%J4KhxSKsjK%)Z&(JGh=~*_VA+PM<xw`ryQ$
ztVTzE2IK}!KlR|ot*fu)RgX4Jdz$<{xaGz*Sy?s980P-in*26Lu3Z<B8xBu6{QTyk
z!sHoD6X#~X^pa{UJs9ko?tl7V-q{NIdN;o2g}d*zMt3g?_}B9!fb-Uq;5}jY#FOh?
zY&IW=VmtWlqa?3&wnD{>%{tX>BHztfkLN0-zIda3K4qSO^5wsi&9$3<d{}6)eBW84
zzZ<=~SEajs4t!pI=~bkkY0$fik0g@pqGlaC{eO<?k;|(;^U8Jq<k-8}JTOelX`Sir
zr8fQy=l@*Se_hTk@cG=!9;E<|?{`k#J;W*Z?agf6zmI1y+`f0We^m&_fp05s70e22
z-u39>vo71DH<#9?&zBHb^ta*CJSmMs@1H$eCCDxq5~wwA^+q<M+Yed}+9`HDQp=7k
zX8FM${q?3ZsG$E1Dd;Cxn}~3T9CP38>8ayer5(AxM^(&ur<C`KWnmsOCtOnO?>KwR
z&uvmth*9+HW5sFJ5=YHCr1X+A+$JCMvCjTfclAl2_VVj_8?A1azk6Tv{NCrk<>$9|
zI4Hcx6yn(Id_n2Y`}wBdPuJy(_3EFQpXwgrbmF|A``g@oPweieY<~D}(IXA*ONJYi
zZ=PF|bdSX?>s;jeXTlaHCMV}u1vDHw_3+1?_H;!ziEE#xpY=`9k5*E*ezI6vZ)ah*
z6z`)WOL7IX`xnpNckrG}mB}l=%Q{Q+y(aW5J}2Fr$Fw@j`Jel#Or5;;kKuXk4!X~l
z)?e^;IusonqW<e}npfGYCi(ot!dqE0Ctc_gzwy7qL{Q)1jms2A{}kK0Tg5*63SQ~W
zI`TC%QcJE$H6ZlA#_NNB0`8t-Tg6mw_Mv!+Q1Qp*Gs^D=wCZM9-Jh)Sb2@vwTp(x4
zEP?emJ?$lqcpk`#yx-Zl@A!^3<C)i=N;5~^ult+xV3NT7O)}S|&AJLR5AJ$*<6h95
zb2(y}s#BWkw|$;=jnyb*cgDI?@wzSSUK6eU#f!GM{q@>q($9bA^ZiRR)~CK;JZZDu
zhE=5RhxDcS_Zt1=ir3H0wS2}I7V*V6@4)Zhxoy4SS|1EoZi|liYwfdJ@1F5NfrN;d
zQ#+^B9+i-&+VG_BHSZQp`<S{ccG;%Tl}+(RsS@kVy3(pm9u!a7J^$4Y`4dm_KU+9V
zcRXI4>$uy|cHIdZqe)*wHM8vBCC_=>Y4|aj_v^%NkrlzoV(e`{3esvr3l68)y=p!F
zBWsV#2bCFlE4KwVEc8s_c^f8OR%NXJe?>Gi&+Ue@j%PShRNL!|uH+pPxNgi}zHWDS
zV)8joo%?A|@@j60PYC3^+CJCxzmn~d3qoy-^|>;)>K5-ft+l$ne4^L#uNP%ACaasz
zk@zKS9DMBlp3724C9ZdhY74nJ=Km{QUc1Hd&lCC2w>Ll4wRrv}e|LnLmr3#If)g93
zEb6zLzNo&ICHbwU%cSEBcV7Ogi@oAtdcx|PRZ&amkCJsuZ#CI)Y05IH+N{?r;7$?w
zmbG!h+Znr}#iE;|mh;Y;duWznS<$jhYxc2y{pt0&B*VG-lr+0pney3`huWTZ7-G&9
z)E0hRBEzB5+PT$$ZS(p3n!Pbb43Fkri5JsNNW3Mpe(k;b@6Jh87fyK>&uhHhVV}FZ
zT-QuYe|y}EmyEx}wVCXLS{$Cg2#<?D>sWmDXnTEy@$cnZb{>rWU)NP6-JAYb&;P*Q
z9FIM_qZzE1yI(Iz4i3EgF)r3sui<vMNd%izOp40#&C_?sx;4eOZa%IgdgSzpd6PCz
zwTPUg#Q!S4=~j2Pou1Leg_afKGO8g5PpjSExb<dlhNyX*MePT{mVDWB6NGnem{zBJ
zBr90b?c8JTcm5g*+t=tOngp}8rIc@e`z6fA_JGU$6V(b_v88V-{|YC(lGp9X({VX?
zPOsSeF`t32wc-62+eDtdeI?T^SI-i6FzwoxF106?EGA(R_KpvK_3c|ytaQcqo587|
z;sT4}KO%+O(wBKhhrC%(w!wenrHN`Q8?<M<{B^E}d+(e!(_LzB7i3sG63NnaFA6vr
ze$?2sqWVJWJ%?aTp4ETe%IXN{1#nr37e^IN&&<2GY;w|um)GRaNQoa@_jj7E<kPbO
ze;(YDT*S<NCgzV>O!7y?pIUFT^8C~Vem~CM*21yz_O<2Ptv1NEu8BLd|5lT{n0nXE
z;F={8`=1=|zn<tH5^}?FlE@Y3DUT1&?sH08wD$VHD!T%ct>!|>sR?KQ)E=LrE?(n$
zkF&JQDdZQ6`?N5v1$!+6&j%|U+M|AD)sNH~9zXm2YbWZl@2=irs~`TYp<ikKzlbAk
z*9?`7zpvps!5va+mLuLgd0WKHhXMWTF9+A8vDomKzwzI8V&Z`llMeb71xjn(DV;X0
zeT!pD%KWcP$#ODF%5qn|);?s`y-3`hN#%gV#Z|UKmHd*&*PrwF8@s@0fr0PyQtekF
zUu`7L?CaLH-?Kv4%|l_o=FuH|>lLmWom{`ojQi3Qb>@RRIro-wr35~A{V@51%N5bh
MX+QWmC72Hb0IoC$CIA2c

literal 0
HcmV?d00001

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 0000000000..b8289d4848
--- /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 <http://www.gnu.org/licenses/>.
+% =========================================================================
+
+% 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 100%
rename from tests/estimation/method_of_moments/AnScho_MoM.mod
rename to tests/estimation/method_of_moments/AnScho/AnScho_MoM.mod
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_MoM_Andreasen.mod b/tests/estimation/method_of_moments/RBC/RBC_MoM_Andreasen.mod
similarity index 100%
rename from tests/estimation/method_of_moments/RBC_MoM_Andreasen.mod
rename to tests/estimation/method_of_moments/RBC/RBC_MoM_Andreasen.mod
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 100%
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
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 100%
rename from tests/estimation/method_of_moments/RBC_MoM_common.inc
rename to tests/estimation/method_of_moments/RBC/RBC_MoM_common.inc
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 100%
rename from tests/estimation/method_of_moments/RBC_MoM_prefilter.mod
rename to tests/estimation/method_of_moments/RBC/RBC_MoM_prefilter.mod
diff --git a/tests/estimation/method_of_moments/RBC_MoM_steady_helper.m b/tests/estimation/method_of_moments/RBC/RBC_MoM_steady_helper.m
similarity index 100%
rename from tests/estimation/method_of_moments/RBC_MoM_steady_helper.m
rename to tests/estimation/method_of_moments/RBC/RBC_MoM_steady_helper.m
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 ba4ef9240b..0000000000
--- 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
-- 
GitLab