diff --git a/meson.build b/meson.build
index bf00091ba6f161a15c84a1db4a7a358444b78568..34c892fb7b7c25b3ed1f0ff9331520ce855ba9c6 100644
--- a/meson.build
+++ b/meson.build
@@ -910,7 +910,47 @@ mod_and_m_tests = [
{ 'test' : [ 'estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod' ],
'extra' : [ 'estimation/method_of_moments/AFVRR/AFVRR_common.inc',
'estimation/method_of_moments/AFVRR/AFVRR_data.mat',
- 'estimation/method_of_moments/AFVRR/AFVRR_steady_helper.m' ] },
+ 'estimation/method_of_moments/AFVRR/AFVRR_steady_helper.m' ] },
+ { 'test' : [ 'estimation/method_of_moments/CET/cet_imh.mod' ],
+ 'extra' : [ 'estimation/method_of_moments/CET/RBC_MoM_common.inc',
+ 'estimation/method_of_moments/CET/cet_matched_irfs_no_interface_workaround.m',
+ 'estimation/method_of_moments/CET/cet_data.mat',
+ 'estimation/method_of_moments/CET/cet_irf_matching_file.m',
+ 'estimation/method_of_moments/CET/cet_model.inc',
+ 'estimation/method_of_moments/CET/cet_original_mode.mat',
+ 'estimation/method_of_moments/CET/cet_steady_helper.m' ] },
+ { 'test' : [ 'estimation/method_of_moments/CET/cet_mle.mod' ],
+ 'extra' : [ 'estimation/method_of_moments/CET/RBC_MoM_common.inc',
+ 'estimation/method_of_moments/CET/cet_matched_irfs_no_interface_workaround.m',
+ 'estimation/method_of_moments/CET/cet_data.mat',
+ 'estimation/method_of_moments/CET/cet_irf_matching_file.m',
+ 'estimation/method_of_moments/CET/cet_model.inc',
+ 'estimation/method_of_moments/CET/cet_original_mode.mat',
+ 'estimation/method_of_moments/CET/cet_steady_helper.m' ] },
+ { 'test' : [ 'estimation/method_of_moments/CET/cet_rwmh.mod' ],
+ 'extra' : [ 'estimation/method_of_moments/CET/RBC_MoM_common.inc',
+ 'estimation/method_of_moments/CET/cet_matched_irfs_no_interface_workaround.m',
+ 'estimation/method_of_moments/CET/cet_data.mat',
+ 'estimation/method_of_moments/CET/cet_irf_matching_file.m',
+ 'estimation/method_of_moments/CET/cet_model.inc',
+ 'estimation/method_of_moments/CET/cet_original_mode.mat',
+ 'estimation/method_of_moments/CET/cet_steady_helper.m' ] },
+ { 'test' : [ 'estimation/method_of_moments/CET/cet_slice.mod' ],
+ 'extra' : [ 'estimation/method_of_moments/CET/RBC_MoM_common.inc',
+ 'estimation/method_of_moments/CET/cet_matched_irfs_no_interface_workaround.m',
+ 'estimation/method_of_moments/CET/cet_data.mat',
+ 'estimation/method_of_moments/CET/cet_irf_matching_file.m',
+ 'estimation/method_of_moments/CET/cet_model.inc',
+ 'estimation/method_of_moments/CET/cet_original_mode.mat',
+ 'estimation/method_of_moments/CET/cet_steady_helper.m' ] },
+ { 'test' : [ 'estimation/method_of_moments/CET/cet_tarb.mod' ],
+ 'extra' : [ 'estimation/method_of_moments/CET/RBC_MoM_common.inc',
+ 'estimation/method_of_moments/CET/cet_matched_irfs_no_interface_workaround.m',
+ 'estimation/method_of_moments/CET/cet_data.mat',
+ 'estimation/method_of_moments/CET/cet_irf_matching_file.m',
+ 'estimation/method_of_moments/CET/cet_model.inc',
+ 'estimation/method_of_moments/CET/cet_original_mode.mat',
+ 'estimation/method_of_moments/CET/cet_steady_helper.m' ] },
{ 'test' : [ 'estimation/conditional-likelihood/1/fs2000_estimation_exact.mod',
'estimation/conditional-likelihood/1/fs2000_estimation_conditional.mod' ],
'extra' : [ 'estimation/fsdat_simul.m' ] },
diff --git a/tests/estimation/method_of_moments/CET/cet_data.mat b/tests/estimation/method_of_moments/CET/cet_data.mat
new file mode 100644
index 0000000000000000000000000000000000000000..5e9f6dd4fa2faa0f3e6cbeac441a29a71e2d2cf2
--- /dev/null
+++ b/tests/estimation/method_of_moments/CET/cet_data.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:1958d890c5c2b01e9b67b0faefeb66e9a0052fcab3a8530d13bb22a8c035f0e2
+size 13553
diff --git a/tests/estimation/method_of_moments/CET/cet_imh.mod b/tests/estimation/method_of_moments/CET/cet_imh.mod
new file mode 100644
index 0000000000000000000000000000000000000000..6aa0c4705a0746e00d24a709c9aa2c1759e9d53c
--- /dev/null
+++ b/tests/estimation/method_of_moments/CET/cet_imh.mod
@@ -0,0 +1,159 @@
+% Functionality testing of Bayesian IRF matching with
+% - independent Metropolis-Hastings
+% - >1 MCMC chains
+% =========================================================================
+% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
+% =========================================================================
+@#include "cet_model.inc"
+
+options_.prior_interval= 0.95;
+
+
+%%%%%%%%%%%%%%%%%%%%%%
+%% NO INTERFACE YET %%
+%%%%%%%%%%%%%%%%%%%%%%
+[M_.matched_irfs, M_.matched_irfs_weights] = cet_matched_irfs_no_interface_workaround(M_.endo_names,M_.exo_names);
+
+
+method_of_moments(mom_method = irf_matching
+%, add_tiny_number_to_cholesky = 1e-14
+%, additional_optimizer_steps = [4]
+%, aim_solver
+%, analytic_jacobian
+%, analytic_standard_errors
+%, bartlett_kernel_lag
+%, bounded_shock_support
+%, brooks_gelman_plotrows
+, cova_compute=0
+%, datafile
+%, dirname=cet_imh_results
+%, dr
+%, dr_cycle_reduction_tol
+%, dr_logarithmic_reduction_maxiter
+%, dr_logarithmic_reduction_tol
+%, drop
+%, first_obs
+%, geweke_interval
+%, graph_format
+%, huge_number
+, irf_matching_file = cet_irf_matching_file
+%, k_order_solver
+%, load_mh_file
+%, load_results_after_load_mh
+%, logdata
+%, lyapunov
+%, lyapunov_complex_threshold
+%, lyapunov_doubling_tol
+%, lyapunov_fixed_point_tol
+, mcmc_jumping_covariance = prior_variance
+%, mh_conf_sig = 0.90
+%, mh_drop = 0.5
+%, mh_init_scale_factor
+%, mh_initialize_from_previous_mcmc
+%, mh_initialize_from_previous_mcmc_directory
+%, mh_initialize_from_previous_mcmc_prior
+%, mh_initialize_from_previous_mcmc_record
+, mh_jscale = 0.5
+, mh_nblocks = 2
+%, mh_posterior_mode_estimation
+%, mh_recover
+, mh_replic=1000
+%, mh_tune_guess = 0.5
+%, mh_tune_jscale = 0.33
+%, mom_burnin
+%, mom_seed
+%, mom_se_tolx
+%, mode_check
+%, mode_check_neighbourhood_size
+%, mode_check_number_of_points
+%, mode_check_symmetric_plots = 0
+, mode_compute = 0
+, mode_file = cet_original_mode
+%, nobs
+%, no_posterior_kernel_density
+%, nodiagnostic
+, nodisplay
+%, nograph
+%, noprint
+%, optim
+%, order
+, penalized_estimator
+, plot_priors = 0
+%, posterior_max_subsample_draws
+%, posterior_sampling_method = 'random_walk_metropolis_hastings'
+%, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+% ,'rand_multivariate_student'
+% ,'student_degrees_of_freedom',4
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% , scale_file
+% )
+%, posterior_sampling_method = 'tailored_random_block_metropolis_hastings'
+%, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+% ,'rand_multivariate_student'
+% ,'student_degrees_of_freedom',4
+% ,'new_block_probability',0.3
+% ,'mode_compute',1
+% ,optim
+% ,'save_tmp_file',1
+% ,scale_file
+% )
+, posterior_sampling_method = 'independent_metropolis_hastings'
+, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+ ,'rand_multivariate_student'
+ ,'student_degrees_of_freedom',4
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% ,scale_file
+ )
+%, posterior_sampling_method = 'slice'
+%, posterior_sampler_options = ('rotated',1
+% ,'mode_files'
+% ,'slice_initialize_with_mode',1
+% ,'initial_step_size',0.8
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% )
+%, prefilter
+%, prior_trunc
+%, pruning
+%, qz_criterium
+%, qz_zero_threshold
+%, raftery_lewis_diagnostics
+%, raftery_lewis_qrs
+%, relative_irf
+%, replic
+%, schur_vec_tol
+%, silent_optimizer
+%, simulation_method = stoch_simul
+%, simulation_multiple
+%, sub_draws
+%, sylvester
+%, sylvester_fixed_point_tol
+%, taper_steps
+%, tex
+%, use_penalized_objective_for_hessian
+%, verbose
+%, weighting_matrix
+%, weighting_matrix_scaling_factor
+%, xls_range
+%, xls_sheet
+);
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/CET/cet_irf_matching_file.m b/tests/estimation/method_of_moments/CET/cet_irf_matching_file.m
new file mode 100644
index 0000000000000000000000000000000000000000..174dc4e44e78aae25dd2de0dfc37893c71c0dbc7
--- /dev/null
+++ b/tests/estimation/method_of_moments/CET/cet_irf_matching_file.m
@@ -0,0 +1,73 @@
+function [modelIrf, check] = cet_irf_matching_file(modelIrf, M_, options_mom_, ys_)
+% Based on replication codes for Christiano, Eichenbaum, Trabandt (2016, Econometrica) - Unemployment and the Business Cycle
+% This showcases how to manipulate model IRFs in a irf_matching_file
+% =========================================================================
+% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
+% =========================================================================
+% - column 1: gdp corresponds to GDPAGG + cumsum(muF)
+% - column 6: real wage corresponds to wAGG + cumsum(muF)
+% - column 7: consumption corresponds to cAGG + cumsum(muF)
+% - column 7: investment corresponds to iAGG + cumsum(muF) + cumsum(mupsiF)
+% - column 8: rel. price investment corresponds to cumsum(pinvestAGG)
+% - column 11: vacancies corresponds to vTotAGG*u
+% - column 12: labor force
+% - column 13: separation rate
+% - column 14: job finding rate corresponds to fAGG*f
+
+% initialize indicator
+check = 0;
+
+modelIrf = 100.*modelIrf;
+
+idmuF = ismember(M_.endo_names,'muF');
+idmupsiF = ismember(M_.endo_names,'mupsiF');
+idGDPAGG = ismember(M_.endo_names,'GDPAGG');
+idpiAGG = ismember(M_.endo_names,'piAGG');
+idRAGG = ismember(M_.endo_names,'RAGG');
+idukAGG = ismember(M_.endo_names,'ukAGG');
+idlAGG = ismember(M_.endo_names,'lAGG');
+idwAGG = ismember(M_.endo_names,'wAGG');
+idcAGG = ismember(M_.endo_names,'cAGG');
+idiAGG = ismember(M_.endo_names,'iAGG');
+idpinvestAGG = ismember(M_.endo_names,'pinvestAGG');
+idpunempAGG = ismember(M_.endo_names,'unempAGG');
+idvTotAGG = ismember(M_.endo_names,'vTotAGG');
+idfAGG = ismember(M_.endo_names,'fAGG');
+
+for jexo=1:M_.exo_nbr
+ if jexo==1
+ mulev_vec = 0;
+ mupsilev_vec = 0;
+ else
+ mulev_vec = cumsum(modelIrf(:,idmuF,jexo));
+ mupsilev_vec = cumsum(modelIrf(:,idmupsiF,jexo));
+ end
+ modelIrf(:,idGDPAGG,jexo) = modelIrf(:,idGDPAGG,jexo) + mulev_vec;
+ modelIrf(:,idpiAGG,jexo) = modelIrf(:,idpiAGG,jexo);
+ modelIrf(:,idRAGG,jexo) = modelIrf(:,idRAGG,jexo);
+ modelIrf(:,idukAGG,jexo) = modelIrf(:,idukAGG,jexo);
+ modelIrf(:,idlAGG,jexo) = modelIrf(:,idlAGG,jexo);
+ modelIrf(:,idwAGG,jexo) = modelIrf(:,idwAGG,jexo) + mulev_vec;
+ modelIrf(:,idcAGG,jexo) = modelIrf(:,idcAGG,jexo) + mulev_vec;
+ modelIrf(:,idiAGG,jexo) = modelIrf(:,idiAGG,jexo) + mulev_vec + mupsilev_vec;
+ modelIrf(:,idpinvestAGG,jexo) = cumsum(modelIrf(:,idpinvestAGG,jexo));
+ modelIrf(:,idpunempAGG,jexo) = modelIrf(:,idpunempAGG,jexo)*M_.params(ismember(M_.param_names,'u'));
+ modelIrf(:,idvTotAGG,jexo) = modelIrf(:,idvTotAGG,jexo);
+ modelIrf(:,idfAGG,jexo) = modelIrf(:,idfAGG,jexo)*M_.params(ismember(M_.param_names,'f'));
+end
+
diff --git a/tests/estimation/method_of_moments/CET/cet_matched_irfs_no_interface_workaround.m b/tests/estimation/method_of_moments/CET/cet_matched_irfs_no_interface_workaround.m
new file mode 100644
index 0000000000000000000000000000000000000000..cd4ad44aca52d9ffee79231fc76c033474ec4502
--- /dev/null
+++ b/tests/estimation/method_of_moments/CET/cet_matched_irfs_no_interface_workaround.m
@@ -0,0 +1,84 @@
+function [matched_irfs, matched_irfs_weights] = cet_matched_irfs_no_interface_workaround(endo_names,exo_names)
+% Based on replication codes for Christiano, Eichenbaum, Trabandt (2016, Econometrica) - Unemployment and the Business Cycle
+% This currently replaces the interface for the IRF Matching capabilities of the method_of_moments toolbox.
+% =========================================================================
+% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
+% =========================================================================
+
+%% settings
+irf_horizon = 15; % horizon of impulse responses to match
+do_monetary_shock_only = 0; % if = 0 all shocks used in estimation
+
+%% load VAR impulse responses from Christiano, Trabandt and Walentin (2010)- Handbook of Monetary Economics Chapter( see main text for reference).
+% IRFFF: impulse responses with respect to monetary policy shock
+% IRFFz: impulse responses with respect to neutral tech shock
+% IRFFu: impulse responses with respect to invest tech shock
+% IRFFFSE, IRFzSE, IRFuSE contain the corresponding estimated standard errors
+% dimensions: rows are periods, columns are variables:
+% - column 1: gdp corresponds to GDPAGG + cumsum(muF)
+% - column 2: inflation corresponds to piAGG
+% - column 3: federal funds rate corresponds to RAGG
+% - column 4: capacity utilization corresponds to ukAGG
+% - column 5: total hours corresponds to lAGG
+% - column 6: real wage corresponds to wAGG + cumsum(muF)
+% - column 7: consumption corresponds to cAGG + cumsum(muF)
+% - column 7: investment corresponds to iAGG + cumsum(muF) + cumsum(mupsiF)
+% - column 8: rel. price investment corresponds to cumsum(pinvestAGG)
+% - column 10: unemployment rate corresponds to unempAGG
+% - column 11: vacancies corresponds to vTotAGG*u
+% - column 12: labor force
+% - column 13: separation rate
+% - column 14: job finding rate corresponds to fAGG*f
+load('cet_data','IRFz','IRFzSE','IRFFF','IRFFFSE','IRFu','IRFuSE');
+
+%% map empirical irf data to a model variable
+% note that any further required transformations or manipulations to model variables (such as cumsum, adding muF and mupsiF)
+% as well as selection of which periods to match occurs in an extra function cet_irf_matching.m
+% if no such function is given then the mapping is exact and the whole horizon will be considered
+
+% irfs with respect to monetary shock
+if do_monetary_shock_only
+ SHOCKNAMES = {'epsR_eps'};
+else
+ SHOCKNAMES = {'epsR_eps', 'muz_eps', 'mupsi_eps'};
+end
+VARNAMES = {'GDPAGG' 'piAGG' 'RAGG' 'ukAGG' 'lAGG' 'wAGG' 'cAGG' 'iAGG' 'pinvestAGG' 'unempAGG' 'vTotAGG' 'fAGG'};
+RESCALE = [1 400 400 1 1 1 1 1 1 100 1 100 ];
+idx=1;
+for jexo = 1:length(SHOCKNAMES)
+ id_shock = strmatch(SHOCKNAMES{jexo},exo_names);
+ if SHOCKNAMES{jexo}=="epsR_eps"
+ IRF = -1*IRFFF(:,[1:11 14]); IRFSE = IRFFFSE(:,[1:11 14]);
+ elseif SHOCKNAMES{jexo}=="muz_eps"
+ IRF = IRFz(:,[1:11 14]); IRFSE = IRFzSE(:,[1:11 14]);
+ elseif SHOCKNAMES{jexo}=="mupsi_eps"
+ IRF = IRFu(:,[1:11 14]); IRFSE = IRFuSE(:,[1:11 14]);
+ end
+ for jvar = 1:length(VARNAMES)
+ id_var = strmatch(VARNAMES{jvar},endo_names);
+ for jirf=1:irf_horizon
+ if IRF(jirf,jvar) ~= 0
+ matched_irfs(idx,:) = {[id_var, id_shock, jirf], IRF(jirf,jvar)/RESCALE(jvar)};
+ matched_irfs_weights(idx,:) = {[id_var, id_shock, jirf], [id_var, id_shock, jirf], 1/(IRFSE(jirf,jvar)/RESCALE(jvar))^2 };
+ idx = idx+1;
+ end
+ end
+ end
+end
+
+end
diff --git a/tests/estimation/method_of_moments/CET/cet_mle.mod b/tests/estimation/method_of_moments/CET/cet_mle.mod
new file mode 100644
index 0000000000000000000000000000000000000000..51c8332f193a497ed86a752d7f6c29415ba33322
--- /dev/null
+++ b/tests/estimation/method_of_moments/CET/cet_mle.mod
@@ -0,0 +1,158 @@
+% Functionality testing of Frequentist IRF matching
+% =========================================================================
+% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
+% =========================================================================
+@#define ML = 1
+@#include "cet_model.inc"
+
+options_.prior_interval= 0.95;
+
+
+%%%%%%%%%%%%%%%%%%%%%%
+%% NO INTERFACE YET %%
+%%%%%%%%%%%%%%%%%%%%%%
+[M_.matched_irfs, M_.matched_irfs_weights] = cet_matched_irfs_no_interface_workaround(M_.endo_names,M_.exo_names);
+
+
+method_of_moments(mom_method = irf_matching
+%, add_tiny_number_to_cholesky = 1e-14
+%, additional_optimizer_steps = [4]
+%, aim_solver
+%, analytic_jacobian
+%, analytic_standard_errors
+%, bartlett_kernel_lag
+%, bounded_shock_support
+%, brooks_gelman_plotrows
+%, cova_compute=0
+%, datafile
+%, dirname=cet_imh_results
+%, dr
+%, dr_cycle_reduction_tol
+%, dr_logarithmic_reduction_maxiter
+%, dr_logarithmic_reduction_tol
+%, drop
+%, first_obs
+%, geweke_interval
+%, graph_format
+%, huge_number
+, irf_matching_file = cet_irf_matching_file
+%, k_order_solver
+%, load_mh_file
+%, load_results_after_load_mh
+%, logdata
+%, lyapunov
+%, lyapunov_complex_threshold
+%, lyapunov_doubling_tol
+%, lyapunov_fixed_point_tol
+%, mcmc_jumping_covariance = prior_variance
+%, mh_conf_sig = 0.90
+%, mh_drop = 0.5
+%, mh_init_scale_factor
+%, mh_initialize_from_previous_mcmc
+%, mh_initialize_from_previous_mcmc_directory
+%, mh_initialize_from_previous_mcmc_prior
+%, mh_initialize_from_previous_mcmc_record
+%, mh_jscale = 0.5
+%, mh_nblocks = 2
+%, mh_posterior_mode_estimation
+%, mh_recover
+%, mh_replic=1000
+%, mh_tune_guess = 0.5
+%, mh_tune_jscale = 0.33
+%, mom_burnin
+%, mom_seed
+%, mom_se_tolx
+%, mode_check
+%, mode_check_neighbourhood_size
+%, mode_check_number_of_points
+%, mode_check_symmetric_plots = 0
+, mode_compute = 4
+%, mode_file = cet_original_mode
+%, nobs
+%, no_posterior_kernel_density
+%, nodiagnostic
+%, nodisplay
+%, nograph
+%, noprint
+%, optim
+%, order
+%, penalized_estimator
+%, plot_priors = 0
+%, posterior_max_subsample_draws
+%, posterior_sampling_method = 'random_walk_metropolis_hastings'
+%, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+% ,'rand_multivariate_student'
+% ,'student_degrees_of_freedom',4
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% , scale_file
+% )
+%, posterior_sampling_method = 'tailored_random_block_metropolis_hastings'
+%, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+% ,'rand_multivariate_student'
+% ,'student_degrees_of_freedom',4
+% ,'new_block_probability',0.3
+% ,'mode_compute',1
+% ,optim
+% ,'save_tmp_file',1
+% ,scale_file
+% )
+%, posterior_sampling_method = 'independent_metropolis_hastings'
+%, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+% ,'rand_multivariate_student'
+% ,'student_degrees_of_freedom',4
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% ,scale_file
+% )
+%, posterior_sampling_method = 'slice'
+%, posterior_sampler_options = ('rotated',1
+% ,'mode_files'
+% ,'slice_initialize_with_mode',1
+% ,'initial_step_size',0.8
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% )
+%, prefilter
+%, prior_trunc
+%, pruning
+%, qz_criterium
+%, qz_zero_threshold
+%, raftery_lewis_diagnostics
+%, raftery_lewis_qrs
+%, relative_irf
+%, replic
+%, schur_vec_tol
+%, silent_optimizer
+, simulation_method = stoch_simul
+%, simulation_multiple
+%, sub_draws
+%, sylvester
+%, sylvester_fixed_point_tol
+%, taper_steps
+%, tex
+%, use_penalized_objective_for_hessian
+%, verbose
+%, weighting_matrix
+%, weighting_matrix_scaling_factor
+%, xls_range
+%, xls_sheet
+);
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/CET/cet_model.inc b/tests/estimation/method_of_moments/CET/cet_model.inc
new file mode 100644
index 0000000000000000000000000000000000000000..c32930f2611fd45e878c432ad620d39a36b34d2d
--- /dev/null
+++ b/tests/estimation/method_of_moments/CET/cet_model.inc
@@ -0,0 +1,751 @@
+/*
+ * Model codes kindly provided by Mathias Trabandt, modified to work with
+ * Dynare 5.0 and newer.
+ *
+ * Please note that the following copyright notice only applies to this Dynare
+ * implementation of the model.
+ */
+
+/*
+ * Copyright © 2023 Dynare Team
+ *
+ * This 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.
+ *
+ * It 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.
+ *
+ * For a copy of the GNU General Public License,
+ * see <https://www.gnu.org/licenses/>.
+ */
+
+//////////////////////////////////////////////////////////////////////////
+//Medium-sized DSGE model of Christiano, Eichenbaum and Trabandt //
+//''Unemployment and Business Cycles'', Econometrica, forthcoming //
+// //
+//Based on parameters: AltOffer(AOB)/Nash bargaining, search/hiring cost//
+//EHL labor market as nested cases. Checkout params.m //
+// //
+//Parameters here are set such that the AOB Model is active. That model //
+//is estimated here and the code reproduces the results. //
+//////////////////////////////////////////////////////////////////////////
+
+
+//Wanted: monetary policy shock is not in period t information set of agents
+//(just as in VAR identified with Choleski decomp.)
+//All other shocks are in period t information set.
+
+//Solution: create two parallel economies which are structually identical. One which is
+//subject to the monetary shock only. This shock is not in period t info set.
+//The other economy has the other shocks with the standard info set.
+
+//variables ending with 'F' are in the economy with full info set
+//variables ending with 'R' are in the restricted info set economy
+//variables ending with 'AGG' are aggregated variables from both economies
+
+//Note that in the 'R' economy, all exogenous variables are set to their steady states, except the monetary shock
+//likewise, in the 'F' economy, the monetary shock is set to zero while all other shocks are active
+
+//Note that both steady states (i.e. 'R' and 'F' economy) are identical
+
+//All Variables in logs.
+
+// Endogenous variables (R economy)
+var wpR wpF wR wF psiR cR RR piR iR pkprimeR FR KR ukR
+kbarR lR xR varthetR phaloR yR JR UR GDPR fR vR
+unempR VR AAR vTotR QR piwR whaloR FwR KwR
+varthetpR;
+
+// Endogenous variables (F economy)
+var psiF cF RF piF iF pkprimeF FF KF ukF
+kbarF lF xF varthetF phaloF yF JF UF GDPF fF
+vF unempF VF AAF vTotF QF piwF whaloF FwF KwF
+muzF mupsiF muF nGF nPHIF nGAMF nDF nKAPF
+varthetpF;
+
+// Endogenous variables (AGG economy)
+var GDPAGG piAGG RAGG ukAGG lAGG wAGG cAGG
+iAGG unempAGG vTotAGG fAGG pinvestAGG;
+
+//Shocks
+varexo epsR_eps muz_eps mupsi_eps;
+
+parameters ydiffdataPercent,idiffdataPercent,alfa,rho,u,Q,sigm,
+betta,lambda,deltak,tau,nuf,etag,kappa,b,xi,D,delta,gamma,Spp
+sigmab sigmaa phi rhoR rpi ry sig_epsR sigmam kappaf DSHARE deltapercent
+recSHAREpercent iota doNash eta doEHL xiw lambdaw AEHL
+f rhomupsi rhomuz sig_mupsi sig_muz thetaG thetaPHI thetaKAP thetaD
+dolagZBAR profy varkappaw pibreve thetaw varkappaf M sigmaL
+tau_SS epsilon_SS upsilon_SS zetac_SS g_SS pibar thetaGAM kappaw
+unemp_SS vTot_SS alp1 alp2 alp3 alp4 bet1 bet2 bet3 s searchSHAREpercent;
+
+model;
+///////////////////////////////////////////////////////////////////////////////////////
+//R economy, monetary policy shock, current realization not in info set////////////////
+///////////////////////////////////////////////////////////////////////////////////////
+//information set in line with Choleski decomposition in VAR
+//set exogenous variables to steady states (except monetary shock, of course)
+#epsilonR=epsilon_SS;
+#upsilonR=upsilon_SS;
+#upsilonR_tp1=upsilon_SS;
+#zetacR=zetac_SS;
+#zetacR_tp1=zetac_SS;
+#gR=g_SS;
+#tauR=tau_SS;
+#mupsiR=STEADY_STATE(mupsiF);
+#mupsiR_tp1=STEADY_STATE(mupsiF);
+#muzR=STEADY_STATE(muzF);
+#muR=alfa/(1-alfa)*STEADY_STATE(mupsiF)+STEADY_STATE(muzF);
+#muR_tp1=STEADY_STATE(muF);
+#nGR=STEADY_STATE(nGF);
+#nPHIR=STEADY_STATE(nPHIF);
+#nGAMR=STEADY_STATE(nGAMF);
+#nDR=STEADY_STATE(nDF);
+#nKAPR=STEADY_STATE(nKAPF);
+
+//conditional expectations based on t-1 info;
+#piR_tp1=EXPECTATION(-1)(piR(+1));
+#FR_tp1=EXPECTATION(-1)(FR(+1));
+#KR_tp1=EXPECTATION(-1)(KR(+1));
+#cR_tp1=EXPECTATION(-1)(cR(+1));
+#psiR_tp1=EXPECTATION(-1)(psiR(+1));
+#pkprimeR_tp1=EXPECTATION(-1)(pkprimeR(+1));
+#iR_tp1=EXPECTATION(-1)(iR(+1));
+#UR_tp1=EXPECTATION(-1)(UR(+1));
+#fR_tp1=EXPECTATION(-1)(fR(+1));
+#wR_tp1=EXPECTATION(-1)(wR(+1));
+#VR_tp1=EXPECTATION(-1)(VR(+1));
+#ukR_tp1=EXPECTATION(-1)(ukR(+1));
+#piwR_tp1=EXPECTATION(-1)(piwR(+1));
+#KwR_tp1=EXPECTATION(-1)(KwR(+1));
+#FwR_tp1=EXPECTATION(-1)(FwR(+1));
+#wpR_tp1=EXPECTATION(-1)(wpR(+1));
+#varthetpR_tp1=EXPECTATION(-1)(varthetpR(+1));
+#AAR_tp1=EXPECTATION(-1)(AAR(+1));
+#RRinfo=EXPECTATION(-1)(RR);
+
+%abbreviations
+#aofukprimeR=sigmab*sigmaa*(exp(ukR))+sigmab*(1-sigmaa);
+#aofukprimeR_tp1=sigmab*sigmaa*(exp(ukR_tp1))+sigmab*(1-sigmaa);
+#aofukR=0.5*sigmab*sigmaa*(exp(ukR))^2+sigmab*(1-sigmaa)*exp(ukR)+sigmab*((sigmaa/2)-1);
+#aofukR_tp1=0.5*sigmab*sigmaa*(exp(ukR_tp1))^2+sigmab*(1-sigmaa)*exp(ukR_tp1)+sigmab*((sigmaa/2)-1);
+#RkR_tp1=log(((exp(ukR_tp1)*aofukprimeR_tp1-aofukR_tp1)+(1-deltak)*exp(pkprimeR_tp1))/(exp(pkprimeR+mupsiR_tp1-piR_tp1)));
+#StildeR=(0.5*(exp(sqrt(Spp)*(exp(iR)/exp(iR(-1))*exp(muR)*exp(mupsiR)-exp(STEADY_STATE(muF))*exp(STEADY_STATE(mupsiF))))+exp(-sqrt(Spp)*(exp(iR)/exp(iR(-1))*exp(muR)*exp(mupsiR)-exp(STEADY_STATE(muF))*exp(STEADY_STATE(mupsiF))))-2));
+#StildeprimeR=(0.5*sqrt(Spp)*(exp(sqrt(Spp)*(exp(iR)/exp(iR(-1))*exp(muR)*exp(mupsiR)-exp(STEADY_STATE(muF))*exp(STEADY_STATE(mupsiF))))-exp(-sqrt(Spp)*(exp(iR)/exp(iR(-1))*exp(muR)*exp(mupsiR)-exp(STEADY_STATE(muF))*exp(STEADY_STATE(mupsiF))))));
+#StildeprimeR_tp1=(0.5*sqrt(Spp)*(exp(sqrt(Spp)*(exp(iR_tp1)/exp(iR)*exp(muR_tp1)*exp(mupsiR_tp1)-exp(STEADY_STATE(muF))*exp(STEADY_STATE(mupsiF))))-exp(-sqrt(Spp)*(exp(iR_tp1)/exp(iR)*exp(muR_tp1)*exp(mupsiR_tp1)-exp(STEADY_STATE(muF))*exp(STEADY_STATE(mupsiF))))));
+#mcR=tauR+alfa*mupsiR+( (1-doEHL)*varthetR+doEHL*wR) +log(nuf*(exp(RRinfo))+1-nuf)-epsilonR-log(1-alfa)-alfa*(kbarR(-1)+ukR-muR-lR-(doEHL*(lambdaw/(lambdaw-1))*whaloR) );
+#pitildewpiR=-piwR+kappaw*piR(-1)+(1-kappaw-varkappaw)*log(pibar)+varkappaw*log(pibreve)+thetaw*STEADY_STATE(muF);
+#pitildewpiR_tp1=-piwR_tp1+kappaw*piR+(1-kappaw-varkappaw)*log(pibar)+varkappaw*log(pibreve)+thetaw*STEADY_STATE(muF);
+#pitildepiR=-piR+kappaf*piR(-1)+(1-kappaf-varkappaf)*log(pibar)+varkappaf*log(pibreve);
+#pitildepiR_tp1=-piR_tp1+kappaf*piR+(1-kappaf-varkappaf)*log(pibar)+varkappaf*log(pibreve);
+
+//R1 - Consumption FOC
+exp(psiR)=exp(zetacR)/(exp(cR)-b*exp(cR(-1))/exp(muR))-betta*b*(exp(zetacR_tp1)/
+(exp(cR_tp1)*exp(muR_tp1)-b*exp(cR)));
+
+//R2 - Bond FOC (Fisher equation)
+exp(psiR)=betta*exp(psiR_tp1)/exp(muR_tp1)*exp(RRinfo)/exp(piR_tp1);
+
+//R3 - Investment FOC
+1=exp(pkprimeR)*exp(upsilonR)*(1-StildeR-StildeprimeR*exp(muR)*exp(mupsiR)*exp(iR)/exp(iR(-1)))+betta*exp(psiR_tp1)/exp(psiR)*exp(pkprimeR_tp1)*exp(upsilonR_tp1)*StildeprimeR_tp1*(exp(iR_tp1)/exp(iR))^2*exp(muR_tp1)*exp(mupsiR_tp1);
+
+//R4 - Capital FOC
+1=betta*exp(psiR_tp1-psiR+RkR_tp1-muR_tp1-piR_tp1);
+
+//R5 - Law of motion for physical capital
+exp(kbarR)=(1-deltak)/(exp(muR)*exp(mupsiR))*exp(kbarR(-1))+exp(upsilonR)*(1-StildeR)*exp(iR);
+
+//R6 - Cost minimization (opt. factor inputs)
+0=aofukprimeR*exp(ukR)*exp(kbarR(-1))/(exp(muR)*exp(mupsiR))-alfa/(1-alfa)*(nuf*exp(RRinfo)+1-nuf)*exp(lR)*((1-doEHL)*(exp(varthetR))+doEHL*( exp(wR)*(exp(whaloR))^(lambdaw/(lambdaw-1))));
+
+//R7 - Production
+exp(yR)=exp(phaloR)^(lambda/(lambda-1))*(exp(epsilonR)*((doEHL*exp(whaloR)^(lambdaw/(lambdaw-1))+(1-doEHL)*1)*exp(lR))^(1-alfa)*
+ (exp(kbarR(-1)+ukR)/(exp(muR)*exp(mupsiR)))^alfa-phi*exp(nPHIR));
+
+//R8 - Resource Constraint
+exp(yR)=exp(gR)*exp(nGR)+exp(cR)+exp(iR)+aofukR*exp(kbarR(-1))/(exp(mupsiR)*exp(muR))
+ +(1-doEHL)*( s*exp(nKAPR)*exp(xR)/exp(QR)*exp(lR(-1)) + kappa*exp(nKAPR)*exp(xR)*exp(lR(-1)) );
+
+//R9 Monetary Policy Rule
+RR=rhoR*RR(-1)+(1-rhoR)*(STEADY_STATE(RR)+rpi*(piR-STEADY_STATE(piF))+ry*(GDPR-STEADY_STATE(GDPF)))-sig_epsR*epsR_eps/400;
+
+//R10 Pricing 1
+exp(FR)-exp(psiR+yR)-betta*xi*exp( pitildepiR_tp1 /(1-lambda)+FR_tp1);
+
+//R11 Pricing 2
+exp(KR)-lambda*exp(psiR+yR+mcR)-betta*xi*exp( pitildepiR_tp1 *lambda/(1-lambda)+KR_tp1);
+
+//R12 Pricing 3
+KR-FR=(1-lambda)*log((1-xi*exp( pitildepiR /(1-lambda)))/(1-xi));
+
+//R13 Price dispersion
+exp(phaloR*lambda/(1-lambda))-(1-xi)^(1-lambda)*(1-xi*exp( pitildepiR /(1-lambda)))^lambda
+ -xi*(exp( pitildepiR +phaloR(-1)))^(lambda/(1-lambda));
+
+//R14 Present value of wages
+doEHL*(exp(wpR)-exp(STEADY_STATE(wpF)))=(1-doEHL)*(-exp(wpR)+exp(wR)+rho*betta*exp(psiR_tp1)/exp(psiR)*exp(wpR_tp1));
+
+//R15 Present value of marginal revenue
+doEHL*(exp(varthetpR)-exp(STEADY_STATE(varthetpF)))=(1-doEHL)*(-exp(varthetpR)+exp(varthetR)+rho*betta*exp(psiR_tp1)/exp(psiR)*exp(varthetpR_tp1));
+
+//R16 Hiring FOC - zero profit/free entry condition
+doEHL*(exp(xR)-exp(STEADY_STATE(xF)))=(1-doEHL)*( -s*exp(nKAPR)/exp(QR) -kappa*exp(nKAPR) + exp(JR) );
+
+//R17 Value of firm
+doEHL*(exp(JR)-exp(STEADY_STATE(JF)))=(1-doEHL)*(-exp(JR)+exp(varthetpR)-exp(wpR));
+
+//R18 Value of work
+doEHL*(exp(VR)-exp(STEADY_STATE(VF)))=(1-doEHL)*(-exp(VR)+exp(wpR)+exp(AAR));
+
+//R19 Present value of worker payoff after separation
+doEHL*(exp(AAR)-exp(STEADY_STATE(AAF)))=(1-doEHL)*(-exp(AAR)+(1-rho)*betta*exp(psiR_tp1)/exp(psiR)*(exp(fR_tp1)*exp(VR_tp1)
++(1-exp(fR_tp1))*exp(UR_tp1))+rho*betta*exp(psiR_tp1)/exp(psiR)*exp(AAR_tp1));
+
+//R20 Unemployment value
+doEHL*(exp(UR)-exp(STEADY_STATE(UF)))=(1-doEHL)*(-exp(UR)+D*exp(nDR)+betta*exp(psiR_tp1)/exp(psiR)*(exp(fR_tp1)*exp(VR_tp1)+(1-exp(fR_tp1))*exp(UR_tp1)));
+
+//R21 Sharing rule
+doEHL*(exp(varthetR)-STEADY_STATE(exp(varthetF)))=(1-doEHL)*(doNash*(exp(VR)-exp(UR)-eta*(exp(JR)+exp(VR)-exp(UR)))+(1-doNash)*(
+
+-exp(JR)+bet1*(exp(VR)-exp(UR))-bet2*exp(nGAMR)*gamma+bet3*(exp(varthetR)-exp(nDR)*D)
+
+));
+
+//R22 GDP
+exp(GDPR)=exp(gR)*exp(nGR)+exp(cR)+exp(iR);
+
+//R23 Unempl. rate
+doEHL*(exp(unempR)-exp(STEADY_STATE(unempF)))=(1-doEHL)*(-exp(unempR)+1-exp(lR));
+
+//R24 Job finding rate
+doEHL*(exp(fR)-exp(STEADY_STATE(fF)))=(1-doEHL)*(-exp(fR)+exp(xR)*exp(lR(-1))/(1-rho*exp(lR(-1))));
+
+//R25 Matching function
+doEHL*(exp(vTotR)-exp(STEADY_STATE(vTotF)))=(1-doEHL)*(-exp(xR)*exp(lR(-1))+sigmam*(1-rho*exp(lR(-1)))^sigm*(exp(vTotR))^(1-sigm));
+
+//R26 Total vacancies
+doEHL*(exp(vR)-exp(STEADY_STATE(vF)))=(1-doEHL)*(-exp(vTotR)+exp(vR)*exp(lR(-1)));
+
+//R27 Vacancy filling rate
+doEHL*(exp(QR)-exp(STEADY_STATE(QF)))=(1-doEHL)*(-exp(QR)+exp(xR)/exp(vR));
+
+//R28 Wage inflation (EHL)
+exp(piwR)=exp(piR)*exp(muR)*exp(wR)/exp(wR(-1));
+
+//R29 Wage dispersion
+doEHL*(exp(whaloR*lambdaw/(1-lambdaw))-(1-xiw)^(1-lambdaw)*(1-xiw*exp(pitildewpiR/(1-lambdaw)))^lambdaw
+ -xiw*(exp(pitildewpiR+whaloR(-1)))^(lambdaw/(1-lambdaw)))=(1-doEHL)*(exp(whaloR)-exp(STEADY_STATE(whaloF)));
+
+//R30 wage setting 1 (EHL)
+doEHL*(-exp(FwR) + exp(psiR)/lambdaw*exp(whaloR)^( lambdaw/(lambdaw-1))*exp(lR) + betta*xiw*(exp(wR_tp1)/exp(wR))*(exp(pitildewpiR_tp1))^( 1/(1-lambdaw) )*exp(FwR_tp1))
+=(1-doEHL)*(exp(FwR)-exp(STEADY_STATE(FwF)));
+
+//R31 Law of motion of employment or EHL wage setting 2
+doEHL*(-exp(KwR) + ( exp(whaloR)^( lambdaw/(lambdaw-1))*exp(lR))^(1+sigmaL) + betta*xiw*(exp(pitildewpiR_tp1))^( lambdaw*(1+sigmaL)/(1-lambdaw) )*exp(KwR_tp1))
+=(1-doEHL)*(-exp(lR)+(rho+exp(xR))*exp(lR(-1)));
+
+//R32 wage setting 3 (EHL)
+doEHL*(1-(1-xiw)*(AEHL*exp(KwR)/(exp(FwR)*exp(wR)))^( 1/( 1-lambdaw*(1+sigmaL) )) - xiw*(exp(pitildewpiR))^(1/(1-lambdaw)))=(1-doEHL)*(exp(KwR)-exp(STEADY_STATE(KwF)));
+
+//////////////////////////////////////////////////////////////////////
+//F economy, other shocks, current realization in info set ///////////
+//////////////////////////////////////////////////////////////////////
+//set not needed exogenous variables to steady states
+#tauF=tau_SS;
+#upsilonF=upsilon_SS;
+#upsilonF_tp1=upsilon_SS;
+#epsilonF=epsilon_SS;
+#zetacF=zetac_SS;
+#zetacF_tp1=zetac_SS;
+#gF=g_SS;
+
+%abbreviations
+#aofukprimeF=sigmab*sigmaa*(exp(ukF))+sigmab*(1-sigmaa);
+#aofukprimeF_tp1=sigmab*sigmaa*(exp(ukF(+1)))+sigmab*(1-sigmaa);
+#aofukF=0.5*sigmab*sigmaa*(exp(ukF))^2+sigmab*(1-sigmaa)*exp(ukF)+sigmab*((sigmaa/2)-1);
+#aofukF_tp1=0.5*sigmab*sigmaa*(exp(ukF(+1)))^2+sigmab*(1-sigmaa)*exp(ukF(+1))+sigmab*((sigmaa/2)-1);
+#RkF_tp1=log(((exp(ukF(+1))*aofukprimeF_tp1-aofukF_tp1)+(1-deltak)*exp(pkprimeF(+1)))/(exp(pkprimeF+mupsiF(+1)-piF(+1))));
+#StildeF=(0.5*(exp(sqrt(Spp)*(exp(iF)/exp(iF(-1))*exp(muF)*exp(mupsiF)-exp(STEADY_STATE(muF))*exp(STEADY_STATE(mupsiF))))+exp(-sqrt(Spp)*(exp(iF)/exp(iF(-1))*exp(muF)*exp(mupsiF)-exp(STEADY_STATE(muF))*exp(STEADY_STATE(mupsiF))))-2));
+#StildeprimeF=(0.5*sqrt(Spp)*(exp(sqrt(Spp)*(exp(iF)/exp(iF(-1))*exp(muF)*exp(mupsiF)-exp(STEADY_STATE(muF))*exp(STEADY_STATE(mupsiF))))-exp(-sqrt(Spp)*(exp(iF)/exp(iF(-1))*exp(muF)*exp(mupsiF)-exp(STEADY_STATE(muF))*exp(STEADY_STATE(mupsiF))))));
+#StildeprimeF_tp1=(0.5*sqrt(Spp)*(exp(sqrt(Spp)*(exp(iF(+1))/exp(iF)*exp(muF(+1))*exp(mupsiF(+1))-exp(STEADY_STATE(muF))*exp(STEADY_STATE(mupsiF))))-exp(-sqrt(Spp)*(exp(iF(+1))/exp(iF)*exp(muF(+1))*exp(mupsiF(+1))-exp(STEADY_STATE(muF))*exp(STEADY_STATE(mupsiF))))));
+#mcF=tauF+alfa*mupsiF+( (1-doEHL)*varthetF+doEHL*wF) +log(nuf*(exp(RF))+1-nuf)-epsilonF-log(1-alfa)-alfa*(kbarF(-1)+ukF-muF-lF-(doEHL*(lambdaw/(lambdaw-1))*whaloF) );
+#pitildewpiF=-piwF+kappaw*piF(-1)+(1-kappaw-varkappaw)*log(pibar)+varkappaw*log(pibreve)+thetaw*STEADY_STATE(muF);
+#pitildewpiF_tp1=-piwF(+1)+kappaw*piF+(1-kappaw-varkappaw)*log(pibar)+varkappaw*log(pibreve)+thetaw*STEADY_STATE(muF);
+#pitildepiF=-piF+kappaf*piF(-1)+(1-kappaf-varkappaf)*log(pibar)+varkappaf*log(pibreve);
+#pitildepiF_tp1=-piF(+1)+kappaf*piF+(1-kappaf-varkappaf)*log(pibar)+varkappaf*log(pibreve);
+
+//F1 - Consumption FOC
+exp(psiF)=exp(zetacF)/(exp(cF)-b*exp(cF(-1))/exp(muF))-betta*b*(exp(zetacF_tp1)/
+(exp(cF(+1))*exp(muF(+1))-b*exp(cF)));
+
+//F2 - Bond FOC (Fisher equation)
+exp(psiF)=betta*exp(psiF(+1))/exp(muF(+1))*exp(RF)/exp(piF(+1));
+
+//F3 - Investment FOC
+1=exp(pkprimeF)*exp(upsilonF)*(1-StildeF-StildeprimeF*exp(muF)*exp(mupsiF)*exp(iF)/exp(iF(-1)))+betta*exp(psiF(+1))/exp(psiF)*exp(pkprimeF(+1))*exp(upsilonF_tp1)*StildeprimeF_tp1*(exp(iF(+1))/exp(iF))^2*exp(muF(+1))*exp(mupsiF(+1));
+
+//F4 - Capital FOC
+1=betta*exp(psiF(+1)-psiF+RkF_tp1-muF(+1)-piF(+1));
+
+//F5 - Law of motion for physical capital
+exp(kbarF)=(1-deltak)/(exp(muF)*exp(mupsiF))*exp(kbarF(-1))+exp(upsilonF)*(1-StildeF)*exp(iF);
+
+//F6 - Cost minimization (opt. factor inputs)
+0=aofukprimeF*exp(ukF)*exp(kbarF(-1))/(exp(muF)*exp(mupsiF))-alfa/(1-alfa)*(nuf*exp(RF)+1-nuf)*exp(lF)*((1-doEHL)*(exp(varthetF))+doEHL*( exp(wF)*(exp(whaloF))^(lambdaw/(lambdaw-1))));
+
+//F7 - Production
+exp(yF)=exp(phaloF)^(lambda/(lambda-1))*(exp(epsilonF)*( (doEHL*exp(whaloF)^(lambdaw/(lambdaw-1))+(1-doEHL)*1)*exp(lF) )^(1-alfa)*
+ (exp(kbarF(-1)+ukF)/(exp(muF)*exp(mupsiF)))^alfa-phi*exp(nPHIF));
+
+//F8 - Resource Constraint
+exp(yF)=exp(gF)*exp(nGF)+exp(cF)+exp(iF)+aofukF*exp(kbarF(-1))/(exp(mupsiF)*exp(muF))
+ +(1-doEHL)*( s*exp(nKAPF)*exp(xF)/exp(QF)*exp(lF(-1)) + kappa*exp(nKAPF)*exp(xF)*exp(lF(-1)) );
+
+
+//F9 Monetary Policy Rule
+RF=rhoR*RF(-1)+(1-rhoR)*(STEADY_STATE(RF)+rpi*(piF-STEADY_STATE(piF))+ry*(GDPF-STEADY_STATE(GDPF)));
+
+//F10 Pricing 1
+exp(FF)-exp(psiF+yF)-betta*xi*exp( pitildepiF_tp1 /(1-lambda)+FF(+1));
+
+//F11 Pricing 2
+exp(KF)-lambda*exp(psiF+yF+mcF)-betta*xi*exp( pitildepiF_tp1*lambda/(1-lambda)+KF(+1));
+
+//F12 Pricing 3
+KF-FF=(1-lambda)*log((1-xi*exp( pitildepiF /(1-lambda)))/(1-xi));
+
+//F13 Price dispersion
+exp(phaloF*lambda/(1-lambda))-(1-xi)^(1-lambda)*(1-xi*exp( pitildepiF /(1-lambda)))^lambda
+ -xi*(exp( pitildepiF +phaloF(-1)))^(lambda/(1-lambda));
+
+//F14 Present value of wages
+doEHL*(exp(wpF)-STEADY_STATE(exp(wpF)))=(1-doEHL)*(-exp(wpF)+exp(wF)+rho*betta*exp(psiF(+1))/exp(psiF)*exp(wpF(+1)));
+
+//F15 Present value of marginal revenue
+doEHL*(exp(varthetpF)-STEADY_STATE(exp(varthetpF)))=(1-doEHL)*(-exp(varthetpF)+exp(varthetF)+rho*betta*exp(psiF(+1))/exp(psiF)*exp(varthetpF(+1)));
+
+//F16 Hiring FOC - zero profit/free entry condition
+doEHL*(exp(xF)-exp(STEADY_STATE(xF)))=(1-doEHL)*( -s*exp(nKAPF)/exp(QF) - kappa*exp(nKAPF) + exp(JF) );
+
+//F17 Firm surplus
+doEHL*(exp(JF)-STEADY_STATE(exp(JF)))=(1-doEHL)*(-exp(JF)+exp(varthetpF)-exp(wpF));
+
+//F18 Value of work
+doEHL*(exp(VF)-exp(STEADY_STATE(VF)))=(1-doEHL)*(-exp(VF)+exp(wpF)+exp(AAF));
+
+//F19 Present value of worker payoff after separation
+doEHL*(exp(AAF)-exp(STEADY_STATE(AAF)))=(1-doEHL)*(-exp(AAF)+(1-rho)*betta*exp(psiF(+1))/exp(psiF)*(exp(fF(+1))*exp(VF(+1))+(1-exp(fF(+1)))*exp(UF(+1)))+rho*betta*exp(psiF(+1))/exp(psiF)*exp(AAF(+1)));
+
+//F20 Unempoyment value
+doEHL*(exp(UF)-exp(STEADY_STATE(UF)))=(1-doEHL)*(-exp(UF)+D*exp(nDF)+betta*exp(psiF(+1))/exp(psiF)*(exp(fF(+1))*exp(VF(+1))+(1-exp(fF(+1)))*exp(UF(+1))));
+
+//R21 Sharing rule
+doEHL*(exp(varthetF)-STEADY_STATE(exp(varthetF)))=(1-doEHL)*(doNash*(exp(VF)-exp(UF)-eta*(exp(JF)+exp(VF)-exp(UF)))+(1-doNash)*(
+
+-exp(JF)+bet1*(exp(VF)-exp(UF))-bet2*exp(nGAMF)*gamma+bet3*(exp(varthetF)-exp(nDF)*D)
+
+));
+
+//F22 GDP
+exp(GDPF)=exp(gF)*exp(nGF)+exp(cF)+exp(iF);
+
+//F23 Unempl. rate
+doEHL*(exp(unempF)-exp(STEADY_STATE(unempF)))=(1-doEHL)*(-exp(unempF)+1-exp(lF));
+
+//F24 Finding rate
+doEHL*(exp(fF)-exp(STEADY_STATE(fF)))=(1-doEHL)*(-exp(fF)+exp(xF)*exp(lF(-1))/(1-rho*exp(lF(-1))));
+
+//F25 Matching function
+doEHL*(exp(vTotF)-exp(STEADY_STATE(vTotF)))=(1-doEHL)*(-exp(xF)*exp(lF(-1))+sigmam*(1-rho*exp(lF(-1)))^sigm*(exp(vTotF))^(1-sigm));
+
+//F26 Total vacancies
+doEHL*(exp(vF)-exp(STEADY_STATE(vF)))=(1-doEHL)*(-exp(vTotF)+exp(vF)*exp(lF(-1)));
+
+//F27 Vacancy filling rate
+doEHL*(exp(QF)-exp(STEADY_STATE(QF)))=(1-doEHL)*(-exp(QF)+exp(xF)/exp(vF));
+
+//F28 Wage inflation
+exp(piwF)=exp(piF)*exp(muF)*exp(wF)/exp(wF(-1));
+
+//F29 Wage dispersion
+doEHL*(exp(whaloF*lambdaw/(1-lambdaw))-(1-xiw)^(1-lambdaw)*(1-xiw*exp(pitildewpiF/(1-lambdaw)))^lambdaw
+ -xiw*(exp(pitildewpiF+whaloF(-1)))^(lambdaw/(1-lambdaw)))=(1-doEHL)*(exp(whaloF)-exp(STEADY_STATE(whaloF)));
+
+//F30 wage setting 1 (EHL)
+doEHL*(-exp(FwF) + exp(psiF)/lambdaw*exp(whaloF)^( lambdaw/(lambdaw-1))*exp(lF) + betta*xiw*(exp(wF(+1))/exp(wF))*(exp(pitildewpiF_tp1))^( 1/(1-lambdaw) )*exp(FwF(+1)) )
+=(1-doEHL)*(exp(FwF)-exp(STEADY_STATE(FwF)));
+
+//F31 Law of motion of employment or EHL wage setting 2
+doEHL*(-exp(KwF) + ( exp(whaloF)^( lambdaw/(lambdaw-1))*exp(lF))^(1+sigmaL) + betta*xiw*( exp(pitildewpiF_tp1) )^( lambdaw*(1+sigmaL)/(1-lambdaw) )*exp(KwF(+1)))
+=(1-doEHL)*(-exp(lF)+(rho+exp(xF))*exp(lF(-1)));
+
+//F32 wage setting 3 (EHL)
+doEHL*(1-(1-xiw)*(AEHL*exp(KwF)/(exp(FwF)*exp(wF)))^( 1/( 1-lambdaw*(1+sigmaL) )) - xiw*( exp(pitildewpiF) )^(1/(1-lambdaw)))=(1-doEHL)*(exp(KwF)-exp(STEADY_STATE(KwF)));
+//////////////////////////
+//Exogenous Variables/////
+//////////////////////////
+
+//E1 Composite technology growth
+muF=alfa/(1-alfa)*mupsiF+muzF;
+
+//E2 Unit root invest. Tech.
+mupsiF=(1-rhomupsi)*STEADY_STATE(mupsiF)+rhomupsi*mupsiF(-1)+sig_mupsi*mupsi_eps/100;
+
+//E3 Unit root neutral Tech.
+muzF=(1-rhomuz)*STEADY_STATE(muzF)+rhomuz*muzF(-1)+sig_muz*muz_eps/100;
+
+//E4 Diffusion of composite technology into gov. spending
+(1-dolagZBAR)*(exp(nGF)-(exp(nGF(-1))/exp(muF))^(1-thetaG))=dolagZBAR*(exp(nGF)-(exp(nGF(-1)))^(1-thetaG)/exp(muF));
+
+//E5 Diffusion of composite technology into fixed costs of production
+nPHIF=nGF;
+
+//E6 Diffusion of composite technology into firm delay costs of bargaining
+nGAMF=nGF;
+
+//E7 Diffusion of composite technology into unemployment benefits
+nDF=nGF;
+
+//E8 Diffusion of composite technology into hiring/search costs
+nKAPF=nGF;
+//////////////////////////////////////////////////////////////////////
+//Aggregating 'F' and 'R' economies, expressed in percent deviations//
+//////////////////////////////////////////////////////////////////////
+//A1
+GDPAGG-STEADY_STATE(GDPF)=GDPR-STEADY_STATE(GDPF)+GDPF-STEADY_STATE(GDPF);
+//A2
+piAGG-STEADY_STATE(piF)=piR-STEADY_STATE(piF)+piF-STEADY_STATE(piF);
+//A3
+RAGG-STEADY_STATE(RF)=RR-STEADY_STATE(RF)+RF-STEADY_STATE(RF);
+//A4
+ukAGG-STEADY_STATE(ukF)=ukR-STEADY_STATE(ukF)+ukF-STEADY_STATE(ukF);
+//A5
+lAGG-STEADY_STATE(lF)=lR-STEADY_STATE(lF)+lF-STEADY_STATE(lF);
+//A6
+wAGG-STEADY_STATE(wF)=wR-STEADY_STATE(wF)+wF-STEADY_STATE(wF);
+//A7
+cAGG-STEADY_STATE(cF)=cR-STEADY_STATE(cF)+cF-STEADY_STATE(cF);
+//A8
+iAGG-STEADY_STATE(iF)=iR-STEADY_STATE(iF)+iF-STEADY_STATE(iF);
+//A9
+unempAGG-STEADY_STATE(unempF)=unempR-STEADY_STATE(unempF)+unempF-STEADY_STATE(unempF);
+//A10
+vTotAGG-STEADY_STATE(vTotF)=vTotR-STEADY_STATE(vTotF)+vTotF-STEADY_STATE(vTotF);
+//A11
+fAGG-STEADY_STATE(fF)=fR-STEADY_STATE(fF)+fF-STEADY_STATE(fF);
+//A12
+pinvestAGG=-mupsiF;
+end;
+/////////////////////////////////
+//End Model //////////
+/////////////////////////////////
+
+
+%read parameters from params.m
+%model switches
+doNash=0; %if =0, alt offer barg. ; if =1, Nash bargaining; note: requires doEHL=0 below
+doEHL=0; %if=0, alt offer barg or Nash barg; if=1, EHL labor market
+@#define do_given_bets=0 // 1: given values for beta's in sharing rule
+
+bet1 = 0.0907; %AOB sharing rule coefficients
+bet2 = 28.9219;
+bet3 = 0.4562;
+
+
+%Labor market parameters
+u=0.055; %unemp. rate
+rho=0.9; %job survival rate
+sigm=0.5570; %matching function share of unemp.
+recSHAREpercent=0.5; %hiring cost para; hiring cost relative to output, in percent
+searchSHAREpercent=0.05; %search cost para; hiring cost relative to output, in percent
+
+DSHARE=0.6682; %unempl. benefits as share of w (replacement ratio)
+Q=0.7; %vacancy filling rate
+
+deltapercent=0.3022; %prob. of barg. session determination
+M=60; %Maximum bargaining rounds per quarter, needs to be an even number!!!!
+
+if M<2, error('M must be at least equal to 2!');end
+if mod(M,2)==1, error('M must be an even number so that the worker makes the last offer!');end
+
+
+%prices
+xi=0.5841; %Calvo prices
+pibar=1.00625; %inflation rate, gross, quarterly
+kappaf=0; %price indexation to past inflation
+varkappaf=1; %price indexation parameter; if kappaf=0 and varkappaf=1 and pibreve=1 -> no indexation at all.
+pibreve=1; %price indexation to value pibreve
+
+lambda=1.4256; %steady state price markup
+nuf=1; %working capital fraction
+tau=1; %steady state markup shifter
+
+%technology and adjustment cost
+alfa=0.2366; %capital share
+deltak=0.025; %depreciation rate of capital
+Spp=13.6684; %second derivative of invest. adjustment costs
+sigmaa=0.0760; %capacity utilization costs
+
+%growth rates of investment and real GDP
+idiffdataPercent=2.9;
+ydiffdataPercent=1.7;
+
+%Preferences
+betta=0.996783170280770; %discount factor households; implies 3% real rate
+b=0.8320; %habit formation in consumption
+
+%Monetary Policy by Taylor rule
+rhoR = 0.8555; %Interest rate smoothing
+rpi = 1.3552; %Coefficient on inflation
+ry = 0.0359; %Coef. on GDP
+
+%Government
+etag=0.2; %Government spending-to-output in ss
+
+%technology diffusion into gov. spending and other parameters of the model
+thetaG =0.0136; %gov spending
+thetaPHI=thetaG; %fixed cost of production
+thetaKAP=thetaG; %hiring cost
+thetaGAM=thetaG; %cost of counteroffer
+thetaD=0.7416; %replacement ratio
+dolagZBAR=1; %if =1, diffusion speed: zbar_t=(zplus_t-1)^thetaj*(zbar_t-1)^(1-thetaj);
+ %if =0, zbar_t=(zplus_t)^thetaj*(zbar_t-1)^(1-thetaj) for j=G,GAM,KAP,BU,PHI
+ %check cet_steadystate.m for further restrictions on
+ %theta parameters!
+
+% steady state profits as share of output
+profy=0;
+
+%EHL parameter
+lambdaw=1.2; %wage markup (EHL)
+xiw=0.75; %Calvo wage stickiness (EHL)
+sigmaL=1; %inv. labor supply elast. (EHL)
+kappaw=0; %wage indexation to past inflation (EHL)
+varkappaw=1; %wage indexation parameter; if kappaw=0 and varkappaw=1 and pibreve=1 and thetaw=0 -> no indexation at all.
+thetaw=0; %wage indexation to technology
+
+
+% standard deviation of shocks (%)
+sig_muz = 0.1362; %unit root neutral tech.
+sig_mupsi = 0.1100; %unit root investment tech.
+sig_epsR = 0.6028; %monetary policy
+
+sig_epsilon = 0; %stationary neutral tech.
+sig_upsilon = 0; %stationary invest. tech.
+sig_g = 0; %gov. spending
+sig_taud = 0; %price markup
+sig_zetac = 0; %cons. preference
+
+%AR(1)
+rhomupsi = 0.7279; %unit root investment tech.
+rhomuz = 0; %unit root neutral tech.
+
+rhoepsilon = 0; %stationary neutral tech.
+rhoupsilon = 0; %stationary invest. tech.
+rhozetac = 0; %cons. preference
+rhog = 0; %gov. spending
+rhotaud = 0; %price markup
+rhosigmam = 0; %matching function shock
+
+
+
+iota=1; %This parameter does not apprear in the model anymore. Just ignore it.
+
+
+shocks;
+var epsR_eps = 1; //monetary policy
+var muz_eps = 1; //neutral tech.
+var mupsi_eps = 1; //invest. tech.
+end;
+
+steady_state_model;
+%% based on cet_steadystate.m, some code is put into cet_steady_helper.m
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%force the below parameters to be identical to thetaG
+thetaPHI=thetaG;%thetaPHI =1; %fixed cost of production
+thetaKAP=thetaG;%thetaKAP =1; %hiring cost
+thetaGAM=thetaG;%thetaGAM =1; %counteroffer cost
+thetaD=thetaG;%thetaD=1; %unemp. benefits
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%first, some easy steady states
+delta=deltapercent/100;
+mu=exp(ydiffdataPercent/100/4);
+mupsi=exp(idiffdataPercent/100/4)/mu;
+muz=mu/mupsi^(alfa/(1-alfa));
+[info,nG,nKAP,nD,nGAM,nPHI] = cet_steady_helper(0,0,0,0,0, dolagZBAR,mu,thetaG,thetaKAP,thetaD,thetaGAM,thetaPHI);
+uk=1;
+aofuk=0;
+Stilde=0;
+Sprimetilde=0;
+l=1-u;
+R=mu*pibar/betta;
+pkprime=1;
+Rk=pibar*mu/betta;
+%Some more complicated steady states
+sigmab=Rk*mupsi*pkprime/pibar-(1-deltak)*pkprime;
+aofukprime=sigmab;
+pitilde=pibar^kappaf*pibar^(1-kappaf-varkappaf)*pibreve^varkappaf;
+mc=1/lambda*(1-betta*xi*(pitilde/pibar)^(lambda/(1-lambda)))/(1-betta*xi*(pitilde/pibar)^(1/(1-lambda))) *((1-xi*(pitilde/pibar)^(1/(1-lambda)))/(1-xi))^(1-lambda);
+phalo=((1-xi*(pitilde/pibar)^(1/(1-lambda)))/(1-xi))^(1-lambda)/((1-xi*(pitilde/pibar)^(lambda/(1-lambda)))/(1-xi))^((1-lambda)/lambda);
+kl=(alfa*(mupsi*mu)^(1-alfa)*mc/sigmab/tau)^(1/(1-alfa));
+%EHL as special case. Obviously, whalo has no meaning otherwise and is set to one!
+piw=mu*pibar;
+pitildew=(pibar)^kappaw*(pibar)^(1-kappaw-varkappaw)*pibreve^varkappaw*mu^thetaw;
+whalo=doEHL*(((1-xiw*(pitildew/piw)^(1/(1-lambdaw)))/(1-xiw))^(1-lambdaw) / ((1-xiw*(pitildew/piw)^(lambdaw/(1-lambdaw)))/(1-xiw))^((1-lambdaw)/lambdaw))+(1-doEHL)*1;
+varthet=(1-alfa)*mc/(tau*(mupsi*mu)^alfa*(nuf*R+1-nuf))*kl^alfa;
+y=mc/( (phalo^(lambda/(1-lambda))-1)*mc+1-profy)*(kl/(mu*mupsi))^alfa*l*whalo^(lambdaw/(lambdaw-1));
+k=kl*l*whalo^(lambdaw/(lambdaw-1));
+phi=((k/(l*whalo^(lambdaw/(lambdaw-1)))/(mupsi*mu))^alfa*l*whalo^(lambdaw/(lambdaw-1))-y*phalo^(lambda/(1-lambda)))/nPHI;
+i=(1-(1-deltak)/(mu*mupsi))*k;
+c=(1-etag)*y-i-(1-doEHL)*recSHAREpercent/100*y-(1-doEHL)*searchSHAREpercent/100*y;
+psi=(c-b*c/mu)^(-1)-betta*b*(c*mu-b*c)^(-1);
+g=etag*y/nG;
+gdp=g*nG+c+i;
+K=lambda*psi*y*mc/(1-betta*xi*(pitilde/pibar)^(lambda/(1-lambda)));
+F=psi*y/(1-betta*xi*(pitilde/pibar)^(1/(1-lambda)));
+%steady states specific to unemployment models, except for real wage (switch)
+x=1-rho;
+f=x*l/(1-rho*l);
+v=x/Q;
+sigmam=x*l*(1-rho*l)^(-sigm)*(l*v)^(sigm-1);
+kappa=recSHAREpercent/100/x/l*y/nKAP;
+s=(searchSHAREpercent/100/((Q^(-1)*x))/l*y)/nKAP;
+J=kappa*nKAP+s*nKAP*Q^(-1);
+varthetp=varthet/(1-rho*betta);
+wp=varthetp-J;
+w=doEHL*(varthet)+(1-doEHL)*(wp*(1-rho*betta)); %EHL switch for wage
+D=DSHARE*w/nD;
+VminusU=(wp-nD*D/(1-rho*betta))/(1-(1-f)*betta*rho)*(1-rho*betta);
+U=(nD*D+betta*f*VminusU)/(1-betta);
+V=VminusU+U;
+A=V-wp;
+%AOB sharing rule
+alp1=1-delta+(1-delta)^M;
+alp2=1-(1-delta)^M;
+alp3=alp2*(1-delta)/delta-alp1;
+alp4=(1-delta)/(2-delta)*alp2/M+1-alp2;
+@#if do_given_bets==0
+ bet1=alp2/alp1;
+ bet2=alp3/alp1;
+ bet3=alp4/alp1;
+@#endif
+gamma=(bet1*(V-U)-J+bet3*(varthet-nD*D))/(nGAM*bet2);
+recruiting=nKAP*kappa*x*l+nKAP*s*(Q^(-1)*x)*l;
+%Worker surplus share under Nash sharing
+eta=(V-U)/(V-U+J);
+%EHL sticky wages
+Kw=((l*whalo^(lambdaw/(lambdaw-1)))^(1+sigmaL))/(1-betta*xiw*(pitildew/piw)^(lambdaw*(1+sigmaL)/(1-lambdaw)));
+Fw=(psi/lambdaw*(l*whalo^(lambdaw/(lambdaw-1))))/(1-betta*xiw*(pitildew/piw)^(1/(1-lambdaw)));
+AEHL=1/(Kw)*(Fw*w)*((1- xiw*( pitildew/piw )^(1/(1-lambdaw)))/(1-xiw))^(( 1-lambdaw*(1+sigmaL) ));
+%parameters for dynare
+tau_SS=log(tau); epsilon_SS=log(1); upsilon_SS=log(1);
+zetac_SS=log(1); g_SS=log(g); unemp_SS=log(1-l); vTot_SS=log(v*l);
+%steady state variables for ys vector (R economy)
+psiR=log(psi); cR=log(c); RR=log(R); piR=log(pibar);
+iR=log(i); pkprimeR=log(pkprime); FR=log(F); KR=log(K);
+ukR=log(uk); kbarR=log(k); lR=log(l); xR=log(x); varthetR=log(varthet);
+phaloR=log(phalo); yR=log(y); VR=log(V); JR=log(J); UR=log(U);
+GDPR=log(gdp); fR=log(f); vR=log(v); unempR=log(u);
+vTotR=log(v*l); QR=log(Q); piwR=log(piw); whaloR=log(whalo);
+FwR=log(Fw); KwR=log(Kw); wR=log(w); wpR=log(wp);
+varthetpR=log(varthetp); AAR=log(A);
+%steady state variables for ys vector (R economy)
+psiF=log(psi); cF=log(c); RF=log(R); piF=log(pibar);
+iF=log(i); pkprimeF=log(pkprime); FF=log(F); KF=log(K);
+ukF=log(uk); kbarF=log(k); lF=log(l); xF=log(x); varthetF=log(varthet);
+phaloF=log(phalo); yF=log(y); VF=log(V); JF=log(J); UF=log(U);
+GDPF=log(gdp); fF=log(f); vF=log(v); unempF=log(u); nDF=log(nD);
+vTotF=log(v*l); QF=log(Q); piwF=log(piw);nPHIF=log(nPHI);wF=log(w);
+whaloF=log(whalo); FwF=log(Fw); KwF=log(Kw); muzF=log(muz);
+mupsiF=log(mupsi); muF=log(mu); nGF=log(nG); nGAMF=log(nGAM); nKAPF=log(nKAP);
+AAF=log(A);wpF=log(wp); varthetpF=log(varthetp);
+%AGG econ
+GDPAGG=log(gdp); piAGG=log(pibar); RAGG=log(R); ukAGG=log(uk); lAGG=log(l);
+wAGG=wF; cAGG=log(c); iAGG=log(i); unempAGG=log(u); vTotAGG=log(v*l);
+fAGG=log(f); pinvestAGG=-log(mupsi);
+
+info = cet_steady_helper(1,Kw,J,V,U); % THIS CHECKS WHETHER STEADY-STATES ARE NONZERO
+end;
+resid;
+steady;
+check;
+
+
+@#ifdef ML
+%% FREQUENTIST ESTIMATION
+estimated_params;
+xi;
+lambda , , 1.001, Inf;
+rhoR , 0.85, -Inf, Inf;
+rpi , 1.45, -Inf, Inf;
+ry , 0.01, -Inf, Inf;
+b , , -Inf, Inf;
+sigmaa , , -Inf, Inf;
+Spp , , -Inf, Inf;
+alfa , , -Inf, 0.6;
+thetaG , 0.10, 0, 1;
+deltapercent , , -Inf, Inf;
+DSHARE , , -Inf, Inf;
+recSHAREpercent , , -Inf, Inf;
+searchSHAREpercent , , -Inf, Inf;
+sigm , , -Inf, Inf;
+sig_epsR , , -Inf, Inf;
+sig_muz , , -Inf, Inf;
+sig_mupsi , , -Inf, Inf;
+rhomupsi , , -Inf, Inf;
+end;
+
+xi = 0.7363;
+lambda = 1.4082;
+rhoR = 0.8463;
+rpi = 1.4566;
+ry = 0.0407;
+b = 0.8196;
+sigmaa = 0.1544;
+Spp = 15.1062;
+alfa = 0.2404;
+thetaG = 0.0502;
+deltapercent = 0.1825;
+DSHARE = 0.4596;
+recSHAREpercent = 0.4849;
+searchSHAREpercent = 0.0654;
+sigm = 0.5403;
+sig_epsR = 0.6163;
+sig_muz = 0.1496;
+sig_mupsi = 0.1150;
+rhomupsi = 0.7314;
+
+estimated_params_init(use_calibration);
+end;
+
+@#else
+
+%% BAYESIAN ESTIMATION
+estimated_params;
+xi , , , , beta_pdf , 0.66, 0.100;
+lambda , , 1.001, , gamma_pdf , 1.20, 0.050;
+rhoR , 0.85, , , beta_pdf , 0.70, 0.150;
+rpi , 1.45, , , gamma_pdf , 1.70, 0.150;
+ry , 0.01, , , gamma_pdf , 0.10, 0.050;
+b , , , , beta_pdf , 0.50, 0.150;
+sigmaa , , , , gamma_pdf , 0.50, 0.300;
+Spp , , , , gamma_pdf , 8.00, 2.000;
+alfa , , , , beta_pdf , 0.33, 0.025;
+thetaG , 0.10, , , beta_pdf , 0.50, 0.200;
+deltapercent , , , , gamma_pdf , 0.50, 0.400;
+DSHARE , , , , beta_pdf , 0.40, 0.100;
+recSHAREpercent , , , , gamma_pdf , 1.00, 0.300;
+searchSHAREpercent , , , , gamma_pdf , 0.10, 0.070;
+sigm , , , , beta_pdf , 0.50, 0.100;
+sig_epsR , , , , gamma_pdf , 0.65, 0.050;
+sig_muz , , , , gamma_pdf , 0.10, 0.050;
+sig_mupsi , , , , gamma_pdf , 0.10, 0.050;
+rhomupsi , , , , beta_pdf , 0.75, 0.100;
+end;
+@#endif
+
+varobs GDPAGG piAGG RAGG ukAGG lAGG wAGG cAGG iAGG pinvestAGG unempAGG vTotAGG fAGG;
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/CET/cet_original_mode.mat b/tests/estimation/method_of_moments/CET/cet_original_mode.mat
new file mode 100644
index 0000000000000000000000000000000000000000..6f804b6dc14568408ebe7d633548f03e808ce0b2
--- /dev/null
+++ b/tests/estimation/method_of_moments/CET/cet_original_mode.mat
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5d6cdcc7cfa32ee9128b0c457798034ffc01024679bcc43f7cbb4ac3ce1840bb
+size 15040
diff --git a/tests/estimation/method_of_moments/CET/cet_rwmh.mod b/tests/estimation/method_of_moments/CET/cet_rwmh.mod
new file mode 100644
index 0000000000000000000000000000000000000000..81076eeb6e25798e396101e6824ad65411774009
--- /dev/null
+++ b/tests/estimation/method_of_moments/CET/cet_rwmh.mod
@@ -0,0 +1,230 @@
+% Functionality testing of Bayesian IRF matching with
+% - Random-Walk Metropolis-Hastings
+% - loading mode files
+% - whether results are stored and can be accessed in different directories
+% - reuse previous MCMC to define covariance of proposal
+% =========================================================================
+% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
+% =========================================================================
+@#include "cet_model.inc"
+
+options_.prior_interval= 0.95;
+
+
+%%%%%%%%%%%%%%%%%%%%%%
+%% NO INTERFACE YET %%
+%%%%%%%%%%%%%%%%%%%%%%
+[M_.matched_irfs, M_.matched_irfs_weights] = cet_matched_irfs_no_interface_workaround(M_.endo_names,M_.exo_names);
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% RUN 1: FIND POSTERIOR MODE USING MODEFILE %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+method_of_moments(mom_method = irf_matching
+%, add_tiny_number_to_cholesky = 1e-14
+, additional_optimizer_steps = [4]
+%, aim_solver
+%, analytic_jacobian
+%, analytic_standard_errors
+%, bartlett_kernel_lag
+%, bounded_shock_support
+%, brooks_gelman_plotrows
+, cova_compute=0
+%, datafile
+, dirname=cet_rwmh_results_1
+%, dr
+%, dr_cycle_reduction_tol
+%, dr_logarithmic_reduction_maxiter
+%, dr_logarithmic_reduction_tol
+%, drop
+%, first_obs
+%, geweke_interval
+%, graph_format
+%, huge_number
+, irf_matching_file = cet_irf_matching_file
+%, k_order_solver
+%, load_mh_file
+%, load_results_after_load_mh
+%, logdata
+%, lyapunov
+%, lyapunov_complex_threshold
+%, lyapunov_doubling_tol
+%, lyapunov_fixed_point_tol
+%, mcmc_jumping_covariance
+, mh_conf_sig = 0.90
+%, mh_drop = 0.5
+%, mh_init_scale_factor
+%, mh_initialize_from_previous_mcmc
+%, mh_initialize_from_previous_mcmc_directory
+%, mh_initialize_from_previous_mcmc_prior
+%, mh_initialize_from_previous_mcmc_record
+%, mh_jscale = 0.5
+%, mh_nblocks = 1
+%, mh_posterior_mode_estimation
+%, mh_recover
+, mh_replic=0
+%, mh_tune_guess = 0.5
+%, mh_tune_jscale = 0.33
+%, mom_burnin
+%, mom_seed
+%, mom_se_tolx
+, mode_check
+%, mode_check_neighbourhood_size
+%, mode_check_number_of_points
+, mode_check_symmetric_plots = 0
+, mode_compute = 1
+, mode_file = cet_original_mode
+%, nobs
+%, no_posterior_kernel_density
+%, nodiagnostic
+%, nodisplay
+%, nograph
+%, noprint
+%, optim
+%, order
+%, penalized_estimator
+, plot_priors = 1
+%, posterior_max_subsample_draws
+%, posterior_sampling_method = 'random_walk_metropolis_hastings'
+%, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+% ,'rand_multivariate_student'
+% ,'student_degrees_of_freedom',4
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% , scale_file
+% )
+%, posterior_sampling_method = 'tailored_random_block_metropolis_hastings'
+%, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+% ,'rand_multivariate_student'
+% ,'student_degrees_of_freedom',4
+% ,'new_block_probability',0.3
+% ,'mode_compute',1
+% ,optim
+% ,'save_tmp_file',1
+% ,scale_file
+% )
+%, posterior_sampling_method = 'independent_metropolis_hastings'
+%, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+% ,'rand_multivariate_student'
+% ,'student_degrees_of_freedom',4
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% ,scale_file
+% )
+%, posterior_sampling_method = 'slice'
+%, posterior_sampler_options = ('rotated',1
+% ,'mode_files'
+% ,'slice_initialize_with_mode',1
+% ,'initial_step_size',0.8
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% )
+%, prefilter
+%, prior_trunc
+%, pruning
+%, qz_criterium
+%, qz_zero_threshold
+%, raftery_lewis_diagnostics
+%, raftery_lewis_qrs
+%, relative_irf
+%, replic
+%, schur_vec_tol
+%, silent_optimizer
+%, simulation_method = stoch_simul
+%, simulation_multiple
+%, sub_draws
+%, sylvester
+%, sylvester_fixed_point_tol
+%, taper_steps
+, tex
+%, use_penalized_objective_for_hessian
+, verbose
+%, weighting_matrix
+%, weighting_matrix_scaling_factor
+%, xls_range
+%, xls_sheet
+);
+close all;
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% RUN 2: LOAD POSTERIOR MODE, COMPUTE HESSIAN AND RUN 1 MCMC CHAIN %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+method_of_moments(mom_method = irf_matching
+, cova_compute=1
+, dirname=cet_rwmh_results_2
+, irf_matching_file = cet_irf_matching_file
+, mh_conf_sig = 0.95
+, mh_drop = 0.5
+, mh_jscale = 0.5
+, mh_nblocks = 1
+, mh_replic=2000
+, mode_compute = 0
+, mode_file = 'cet_rwmh_results_1/method_of_moments/cet_rwmh_mode'
+, plot_priors = 0
+, posterior_sampling_method = 'random_walk_metropolis_hastings'
+, posterior_sampler_options = ('proposal_distribution'
+ ,'rand_multivariate_normal'
+% ,'rand_multivariate_student'
+% ,'student_degrees_of_freedom',4
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% , scale_file
+ )
+, tex
+);
+close all;
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% RUN 3: DEFINE COVARIANCE OF PROPOSAL FROM PREVIOUS MCMC %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+method_of_moments(mom_method = irf_matching
+, dirname=cet_rwmh_results_2
+, irf_matching_file = cet_irf_matching_file
+, load_mh_file
+, mh_conf_sig = 0.95
+, mh_drop = 0.5
+, mh_jscale = 0.5
+, mh_nblocks = 1
+, mh_replic=600
+, mode_compute = 0
+, plot_priors = 0
+, posterior_sampling_method = 'random_walk_metropolis_hastings'
+, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+ ,'rand_multivariate_student'
+ ,'student_degrees_of_freedom',4
+ ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% , scale_file
+ )
+, tex
+);
+
+
+%%%%%%%%%%%
+%% LATEX %%
+%%%%%%%%%%%
+collect_latex_files;
+if system(['pdflatex -halt-on-error -interaction=batchmode ' M_.fname '_TeX_binder.tex'])
+ error('TeX-File did not compile.')
+end
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/CET/cet_slice.mod b/tests/estimation/method_of_moments/CET/cet_slice.mod
new file mode 100644
index 0000000000000000000000000000000000000000..a03d4b90579a8e0ab9e5d793f12e2dc37628f121
--- /dev/null
+++ b/tests/estimation/method_of_moments/CET/cet_slice.mod
@@ -0,0 +1,206 @@
+% Functionality testing of Bayesian IRF matching with
+% - slice
+% - rotated slice with use_mh_covariance_matrix
+% - rotated slice with slice_initialize_with_mode
+% - reuse previous MCMC
+% =========================================================================
+% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
+% =========================================================================
+@#include "cet_model.inc"
+
+options_.prior_interval= 0.95;
+
+
+%%%%%%%%%%%%%%%%%%%%%%
+%% NO INTERFACE YET %%
+%%%%%%%%%%%%%%%%%%%%%%
+[M_.matched_irfs, M_.matched_irfs_weights] = cet_matched_irfs_no_interface_workaround(M_.endo_names,M_.exo_names);
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% RUN 1: DEFAULT SLICE %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%
+method_of_moments(mom_method = irf_matching
+%, add_tiny_number_to_cholesky = 1e-14
+%, additional_optimizer_steps = [4]
+%, aim_solver
+%, analytic_jacobian
+%, analytic_standard_errors
+%, bartlett_kernel_lag
+%, bounded_shock_support
+%, brooks_gelman_plotrows
+%, cova_compute=0
+%, datafile
+%, dirname=cet_slice_default
+%, dr
+%, dr_cycle_reduction_tol
+%, dr_logarithmic_reduction_maxiter
+%, dr_logarithmic_reduction_tol
+%, drop
+%, first_obs
+%, geweke_interval
+%, graph_format
+%, huge_number
+, irf_matching_file = cet_irf_matching_file
+%, k_order_solver
+%, load_mh_file
+%, load_results_after_load_mh
+%, logdata
+%, lyapunov
+%, lyapunov_complex_threshold
+%, lyapunov_doubling_tol
+%, lyapunov_fixed_point_tol
+%, mcmc_jumping_covariance
+, mh_conf_sig = 0.90
+%, mh_drop = 0.5
+%, mh_init_scale_factor
+%, mh_initialize_from_previous_mcmc
+%, mh_initialize_from_previous_mcmc_directory
+%, mh_initialize_from_previous_mcmc_prior
+%, mh_initialize_from_previous_mcmc_record
+%, mh_jscale = 0.5
+, mh_nblocks = 1
+%, mh_posterior_mode_estimation
+%, mh_recover
+, mh_replic=80
+%, mh_tune_guess = 0.5
+%, mh_tune_jscale = 0.33
+%, mom_burnin
+%, mom_seed
+%, mom_se_tolx
+%, mode_check
+%, mode_check_neighbourhood_size
+%, mode_check_number_of_points
+%, mode_check_symmetric_plots = 0
+%, mode_compute = 1
+%, mode_file = cet_original_mode
+%, nobs
+%, no_posterior_kernel_density
+%, nodiagnostic
+%, nodisplay
+%, nograph
+%, noprint
+%, optim
+%, order
+%, penalized_estimator
+, plot_priors = 0
+%, posterior_max_subsample_draws
+%, posterior_sampling_method = 'random_walk_metropolis_hastings'
+%, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+% ,'rand_multivariate_student'
+% ,'student_degrees_of_freedom',4
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% , scale_file
+% )
+%, posterior_sampling_method = 'tailored_random_block_metropolis_hastings'
+%, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+% ,'rand_multivariate_student'
+% ,'student_degrees_of_freedom',4
+% ,'new_block_probability',0.3
+% ,'mode_compute',1
+% ,optim
+% ,'save_tmp_file',1
+% ,scale_file
+% )
+%, posterior_sampling_method = 'independent_metropolis_hastings'
+%, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+% ,'rand_multivariate_student'
+% ,'student_degrees_of_freedom',4
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% ,scale_file
+% )
+, posterior_sampling_method = 'slice'
+%, posterior_sampler_options = ('rotated',1
+% ,'mode_files'
+% ,'slice_initialize_with_mode',1
+% ,'initial_step_size',0.8
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',0
+% )
+%, prefilter
+%, prior_trunc
+%, pruning
+%, qz_criterium
+%, qz_zero_threshold
+%, raftery_lewis_diagnostics
+%, raftery_lewis_qrs
+%, relative_irf
+%, replic
+%, schur_vec_tol
+%, silent_optimizer
+, simulation_method = stoch_simul
+%, simulation_multiple
+%, sub_draws
+%, sylvester
+%, sylvester_fixed_point_tol
+%, taper_steps
+%, tex
+%, use_penalized_objective_for_hessian
+%, verbose
+%, weighting_matrix
+%, weighting_matrix_scaling_factor
+%, xls_range
+%, xls_sheet
+);
+close all;
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% RUN 2: ROTATED SLICE WITH USE_MH_COVARIANCE_MATRIX %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+method_of_moments(mom_method = irf_matching
+, irf_matching_file = cet_irf_matching_file
+, mh_nblocks = 1
+, mh_replic=20
+, plot_priors = 0
+, load_mh_file
+, posterior_sampling_method = 'slice'
+, posterior_sampler_options = ('rotated',1
+% ,'mode_files'
+% ,'slice_initialize_with_mode',1
+% ,'initial_step_size',0.8
+ ,'use_mh_covariance_matrix',1
+ ,'save_tmp_file',0
+ )
+);
+close all;
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% RUN 3: ROTATED SLICE WITH SLICE_INITIALIZE_WITH_MODE %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+method_of_moments(mom_method = irf_matching
+, irf_matching_file = cet_irf_matching_file
+, mh_nblocks = 1
+, mh_replic=20
+, plot_priors = 0
+, mode_compute = 1
+, posterior_sampling_method = 'slice'
+, posterior_sampler_options = ('rotated',1
+% ,'mode_files'
+ ,'slice_initialize_with_mode',1
+% ,'initial_step_size',0.8
+% ,'use_mh_covariance_matrix',1
+ ,'save_tmp_file',0
+ )
+);
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/CET/cet_steady_helper.m b/tests/estimation/method_of_moments/CET/cet_steady_helper.m
new file mode 100644
index 0000000000000000000000000000000000000000..e137d83d32934b0c0138be43080b3135be959cb7
--- /dev/null
+++ b/tests/estimation/method_of_moments/CET/cet_steady_helper.m
@@ -0,0 +1,50 @@
+function [info,nG,nKAP,nD,nGAM,nPHI] = cet_steady_helper(justCheck,Kw,J,V,U, dolagZBAR,mu,thetaG,thetaKAP,thetaD,thetaGAM,thetaPHI)
+% Based on cet_steadystate.m of the replication codes for
+% Christiano, Eichenbaum, Trabandt (2016, Econometrica) - Unemployment and the Business Cycle;
+% slightly modified such that most code is now in a steady_state_model block
+% and only the two if statements are computed in this helper function.
+% =========================================================================
+% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
+% =========================================================================
+
+info = 0;
+if justCheck
+ if (Kw<0)
+ disp('could not compute steadystate: Kw<0');
+ elseif (J<0)
+ disp('could not compute steadystate: J<0');
+ elseif (V-U<0)
+ disp('could not compute steadystate: V-U<0')
+ info=1;
+ end
+ return
+else
+ if dolagZBAR==1
+ nG=mu^(-1/thetaG);
+ nKAP=mu^(-1/thetaKAP);
+ nD=mu^(-1/thetaD);
+ nGAM=mu^(-1/thetaGAM);
+ nPHI=mu^(-1/thetaPHI);
+ else
+ nG=(1/mu)^((1-thetaG)/thetaG);
+ nKAP=(1/mu)^((1-thetaKAP)/thetaKAP);
+ nD=(1/mu)^((1-thetaD)/thetaD);
+ nGAM=(1/mu)^((1-thetaGAM)/thetaGAM);
+ nPHI=(1/mu)^((1-thetaPHI)/thetaPHI);
+ end
+end
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/CET/cet_tarb.mod b/tests/estimation/method_of_moments/CET/cet_tarb.mod
new file mode 100644
index 0000000000000000000000000000000000000000..63f49b4426adc510d107a596a298980394f141fe
--- /dev/null
+++ b/tests/estimation/method_of_moments/CET/cet_tarb.mod
@@ -0,0 +1,181 @@
+% Functionality testing of Bayesian IRF matching with
+% - Tailored-Randomized-Block Metropolis-Hastings
+% - reuse previous MCMC mode for optimization
+% =========================================================================
+% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
+% =========================================================================
+@#include "cet_model.inc"
+
+options_.prior_interval= 0.95;
+
+
+%%%%%%%%%%%%%%%%%%%%%%
+%% NO INTERFACE YET %%
+%%%%%%%%%%%%%%%%%%%%%%
+[M_.matched_irfs, M_.matched_irfs_weights] = cet_matched_irfs_no_interface_workaround(M_.endo_names,M_.exo_names);
+
+
+method_of_moments(mom_method = irf_matching
+%, add_tiny_number_to_cholesky = 1e-14
+, additional_optimizer_steps = [4]
+%, aim_solver
+%, analytic_jacobian
+%, analytic_standard_errors
+%, bartlett_kernel_lag
+%, bounded_shock_support
+%, brooks_gelman_plotrows
+, cova_compute=1
+%, datafile
+, dirname=cet_tarb_results
+%, dr
+%, dr_cycle_reduction_tol
+%, dr_logarithmic_reduction_maxiter
+%, dr_logarithmic_reduction_tol
+%, drop
+%, first_obs
+%, geweke_interval
+%, graph_format
+%, huge_number
+, irf_matching_file = cet_irf_matching_file
+%, k_order_solver
+%, load_mh_file
+%, load_results_after_load_mh
+%, logdata
+%, lyapunov
+%, lyapunov_complex_threshold
+%, lyapunov_doubling_tol
+%, lyapunov_fixed_point_tol
+%, mcmc_jumping_covariance
+, mh_conf_sig = 0.90
+%, mh_drop = 0.5
+%, mh_init_scale_factor
+%, mh_initialize_from_previous_mcmc
+%, mh_initialize_from_previous_mcmc_directory
+%, mh_initialize_from_previous_mcmc_prior
+%, mh_initialize_from_previous_mcmc_record
+, mh_jscale = 0.5
+, mh_nblocks = 1
+%, mh_posterior_mode_estimation
+%, mh_recover
+, mh_replic=25
+%, mh_tune_guess = 0.5
+%, mh_tune_jscale = 0.33
+%, mom_burnin
+%, mom_seed
+%, mom_se_tolx
+, mode_check
+%, mode_check_neighbourhood_size
+%, mode_check_number_of_points
+%, mode_check_symmetric_plots = 0
+, mode_compute = 1
+, mode_file = cet_original_mode
+%, nobs
+%, no_posterior_kernel_density
+%, nodiagnostic
+%, nodisplay
+%, nograph
+%, noprint
+%, optim
+%, order
+%, penalized_estimator
+, plot_priors = 0
+%, posterior_max_subsample_draws
+%, posterior_sampling_method = 'random_walk_metropolis_hastings'
+%, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+% ,'rand_multivariate_student'
+% ,'student_degrees_of_freedom',4
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% , scale_file
+% )
+, posterior_sampling_method = 'tailored_random_block_metropolis_hastings'
+, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+ ,'rand_multivariate_student'
+ ,'student_degrees_of_freedom',4
+ ,'new_block_probability',0.25
+ ,'mode_compute',4
+% ,optim
+% ,'save_tmp_file',1
+% ,scale_file
+ )
+%, posterior_sampling_method = 'independent_metropolis_hastings'
+%, posterior_sampler_options = ('proposal_distribution'
+% ,'rand_multivariate_normal'
+% ,'rand_multivariate_student'
+% ,'student_degrees_of_freedom',4
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% ,scale_file
+% )
+%, posterior_sampling_method = 'slice'
+%, posterior_sampler_options = ('rotated',1
+% ,'mode_files'
+% ,'slice_initialize_with_mode',1
+% ,'initial_step_size',0.8
+% ,'use_mh_covariance_matrix',1
+% ,'save_tmp_file',1
+% )
+%, prefilter
+%, prior_trunc
+%, pruning
+%, qz_criterium
+%, qz_zero_threshold
+%, raftery_lewis_diagnostics
+%, raftery_lewis_qrs
+%, relative_irf
+%, replic
+%, schur_vec_tol
+%, silent_optimizer
+%, simulation_method = stoch_simul
+%, simulation_multiple
+%, sub_draws
+%, sylvester
+%, sylvester_fixed_point_tol
+%, taper_steps
+, tex
+%, use_penalized_objective_for_hessian
+%, verbose
+%, weighting_matrix
+%, weighting_matrix_scaling_factor
+%, xls_range
+%, xls_sheet
+);
+
+
+method_of_moments(mom_method = irf_matching
+, additional_optimizer_steps = [1]
+, cova_compute=1
+, dirname=cet_tarb_results
+, irf_matching_file = cet_irf_matching_file
+, mh_conf_sig = 0.90
+, mh_replic=0
+, mode_check
+, mode_compute = 4
+, mode_file = 'cet_tarb_results/method_of_moments/cet_tarb_mh_mode'
+, plot_priors = 0
+);
+
+%%%%%%%%%%%
+%% LATEX %%
+%%%%%%%%%%%
+collect_latex_files;
+if system(['pdflatex -halt-on-error -interaction=batchmode ' M_.fname '_TeX_binder.tex'])
+ error('TeX-File did not compile.')
+end
\ No newline at end of file