diff --git a/matlab/+mom/default_option_mom_values.m b/matlab/+mom/default_option_mom_values.m
index 82d05b86da9a4423f0fdfb911b79abf82542d90c..90c383d291d0d5fb07bcc1ec91cf104c9e892bed 100644
--- a/matlab/+mom/default_option_mom_values.m
+++ b/matlab/+mom/default_option_mom_values.m
@@ -1,5 +1,5 @@
-function options_mom_ = default_option_mom_values(options_mom_, options_, dname, doBayesianEstimation)
-% function options_mom_ = default_option_mom_values(options_mom_, options_, dname, doBayesianEstimation)
+function options_mom_ = default_option_mom_values(options_mom_, options_, fname, doBayesianEstimation)
+% function options_mom_ = default_option_mom_values(options_mom_, options_, fname, doBayesianEstimation)
 % -------------------------------------------------------------------------
 % Returns structure containing the options for method_of_moments command.
 % Note that options_mom_ is local and contains default and user-specified
@@ -15,7 +15,7 @@ function options_mom_ = default_option_mom_values(options_mom_, options_, dname,
 % INPUTS
 %  o options_mom_:         [structure]  all user-specified settings (from the method_of_moments command)
 %  o options_:             [structure]  Dynare global options
-%  o dname:                [string]     default name of directory to store results (usually name of mod file)
+%  o fname:                [string]     name of mod file (used as default for dirname)
 %  o doBayesianEstimation  [boolean]    indicator whether we do Bayesian estimation
 % -------------------------------------------------------------------------
 % OUTPUTS
@@ -82,7 +82,7 @@ options_mom_.bayesian_irf = false;
 % -------------------------------------------------------------------------
 
 % common settings
-options_mom_ = set_default_option(options_mom_,'dirname',dname);         % specify directory in which to store estimation output
+options_mom_ = set_default_option(options_mom_,'dirname',fname);         % specify directory in which to store estimation output
 options_mom_ = set_default_option(options_mom_,'graph_format','eps');    % specify the file format(s) for graphs saved to disk
 options_mom_ = set_default_option(options_mom_,'nodisplay',false);       % do not display the graphs, but still save them to disk
 options_mom_ = set_default_option(options_mom_,'nograph',false);         % do not create graphs (which implies that they are not saved to the disk nor displayed)
@@ -95,6 +95,7 @@ if doBayesianEstimation
 end
 
 % specific method_of_moments settings
+options_mom_.mom = set_default_option(options_mom_.mom,'simulation_method','STOCH_SIMUL'); % simulation method used to compute moments or IRFs
 if strcmp(mom_method,'GMM') || strcmp(mom_method,'SMM')
     options_mom_.mom = set_default_option(options_mom_.mom,'bartlett_kernel_lag',20);                   % bandwith in optimal weighting matrix
     options_mom_.mom = set_default_option(options_mom_.mom,'penalized_estimator',false);                % include deviation from prior mean as additional moment restriction and use prior precision as weights
@@ -113,14 +114,15 @@ if strcmp(mom_method,'GMM')
     options_mom_.mom = set_default_option(options_mom_.mom,'analytic_standard_errors',false); % compute standard errors numerically (0) or analytically (1). Analytical derivatives are only available for GMM.
 end
 if strcmp(mom_method,'IRF_MATCHING')
-    options_mom_.mom = set_default_option(options_mom_.mom,'simulation_method','STOCH_SIMUL'); % simulation method used to compute IRFs
     if ~isfield(options_mom_.mom,'irf_matching_file') 
         options_mom_.mom.irf_matching_file = [];  % irf_matching file enables to transform model IRFs before matching them to data IRFs
     end
     options_mom_.mom.irf_matching_file = set_default_option(options_mom_.mom.irf_matching_file,'name','');    
-    options_mom_ = set_default_option(options_mom_,'add_tiny_number_to_cholesky',1e-14); % add tiny number to Cholesky factor to avoid numerical problems when computing IRFs
-    options_mom_ = set_default_option(options_mom_,'drop',100); % truncation / burnin for order>1 irf simulations
-    options_mom_ = set_default_option(options_mom_,'relative_irf',false); % requests the computation of normalized IRFs    
+    if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+        options_mom_ = set_default_option(options_mom_,'add_tiny_number_to_cholesky',1e-14); % add tiny number to Cholesky factor to avoid numerical problems when computing IRFs    
+        options_mom_ = set_default_option(options_mom_,'drop',100); % truncation / burnin for order>1 irf simulations
+        options_mom_ = set_default_option(options_mom_,'relative_irf',false); % requests the computation of normalized IRFs
+    end
 end
 
 % data related options
@@ -155,37 +157,54 @@ options_mom_ = set_default_option(options_mom_,'mode_file','');
 options_mom_ = set_default_option(options_mom_,'cova_compute',true);               % 1: computed covariance via Hessian after the computation of the mode, 0: turn off computation of covariance matrix
 
 % perturbation related
-options_mom_ = set_default_option(options_mom_,'order',1);                              % order of Taylor approximation in perturbation
-if strcmp(mom_method,'IRF_MATCHING')                                                    % number of simulated series used to compute IRFs
-    if options_mom_.order == 1
-        options_mom_ = set_default_option(options_mom_,'replic',1);
-    else
-        options_mom_ = set_default_option(options_mom_,'replic',50);
+if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+    options_mom_ = set_default_option(options_mom_,'order',1);                              % order of Taylor approximation in perturbation
+    if strcmp(mom_method,'IRF_MATCHING')                                                    % number of simulated series used to compute IRFs
+        if options_mom_.order == 1
+            options_mom_ = set_default_option(options_mom_,'replic',1);
+        else
+            options_mom_ = set_default_option(options_mom_,'replic',50);
+        end
     end
+    options_mom_ = set_default_option(options_mom_,'pruning',false);                        % use pruned state space system at order>1
+    options_mom_ = set_default_option(options_mom_,'aim_solver',false);                     % use AIM algorithm to compute perturbation approximation instead of mjdgges
+    options_mom_ = set_default_option(options_mom_,'k_order_solver',false);                 % use k_order_perturbation instead of mjdgges
+    options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction',false);             % use cycle reduction algorithm to solve the polynomial equation for retrieving the coefficients associated to the endogenous variables in the decision rule
+    options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction_tol',1e-7);          % convergence criterion used in the cycle reduction algorithm
+    options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction',false);       % use logarithmic reduction algorithm to solve the polynomial equation for retrieving the coefficients associated to the endogenous variables in the decision rule
+    options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_maxiter',100); % maximum number of iterations used in the logarithmic reduction algorithm
+    options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_tol',1e-12);   % convergence criterion used in the cycle reduction algorithm
+    options_mom_ = set_default_option(options_mom_,'qz_criterium',1-1e-6);                  % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
+                                                                                            % if there are no unit roots one can use 1.0 (or slightly below) which we set as default; if they are possible, you may have have multiple unit roots and the accuracy decreases when computing the eigenvalues in lyapunov_symm
+                                                                                            % Note that unit roots are only possible at first-order, at higher order we set it to 1 in pruned_state_space_system and focus only on stationary observables.
+    options_mom_ = set_default_option(options_mom_,'qz_zero_threshold',1e-6);               % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
+    options_mom_ = set_default_option(options_mom_,'schur_vec_tol',1e-11);                  % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix.
 end
-options_mom_ = set_default_option(options_mom_,'pruning',false);                        % use pruned state space system at order>1
-options_mom_ = set_default_option(options_mom_,'aim_solver',false);                     % use AIM algorithm to compute perturbation approximation instead of mjdgges
-options_mom_ = set_default_option(options_mom_,'k_order_solver',false);                 % use k_order_perturbation instead of mjdgges
-options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction',false);             % use cycle reduction algorithm to solve the polynomial equation for retrieving the coefficients associated to the endogenous variables in the decision rule
-options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction_tol',1e-7);          % convergence criterion used in the cycle reduction algorithm
-options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction',false);       % use logarithmic reduction algorithm to solve the polynomial equation for retrieving the coefficients associated to the endogenous variables in the decision rule
-options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_maxiter',100); % maximum number of iterations used in the logarithmic reduction algorithm
-options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_tol',1e-12);   % convergence criterion used in the cycle reduction algorithm
-options_mom_ = set_default_option(options_mom_,'qz_criterium',1-1e-6);                  % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
-                                                                                        % if there are no unit roots one can use 1.0 (or slightly below) which we set as default; if they are possible, you may have have multiple unit roots and the accuracy decreases when computing the eigenvalues in lyapunov_symm
-                                                                                        % Note that unit roots are only possible at first-order, at higher order we set it to 1 in pruned_state_space_system and focus only on stationary observables.
-options_mom_ = set_default_option(options_mom_,'qz_zero_threshold',1e-6);               % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
-options_mom_ = set_default_option(options_mom_,'schur_vec_tol',1e-11);                  % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix.
 
-% numerical algorithms
-options_mom_ = set_default_option(options_mom_,'lyapunov_db',false);                    % doubling algorithm (disclyap_fast) to solve Lyapunov equation to compute variance-covariance matrix of state variables
-options_mom_ = set_default_option(options_mom_,'lyapunov_fp',false);                    % fixed-point algorithm to solve Lyapunov equation to compute variance-covariance matrix of state variables
-options_mom_ = set_default_option(options_mom_,'lyapunov_srs',false);                   % square-root-solver (dlyapchol) algorithm to solve Lyapunov equation to compute variance-covariance matrix of state variables
-options_mom_ = set_default_option(options_mom_,'lyapunov_complex_threshold',1e-15);     % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
-options_mom_ = set_default_option(options_mom_,'lyapunov_fixed_point_tol',1e-10);       % convergence criterion used in the fixed point Lyapunov solver
-options_mom_ = set_default_option(options_mom_,'lyapunov_doubling_tol',1e-16);          % convergence criterion used in the doubling algorithm
-options_mom_ = set_default_option(options_mom_,'sylvester_fp',false);                   % determines whether to use fixed point algorihtm to solve Sylvester equation (gensylv_fp), faster for large scale models
-options_mom_ = set_default_option(options_mom_,'sylvester_fixed_point_tol',1e-12);      % convergence criterion used in the fixed point Sylvester solver
+% perfect foresight related
+if strcmp(options_mom_.mom.simulation_method,'PERFECT_FORESIGHT')
+    options_mom_ = set_default_option(options_mom_,'periods',options_.periods);
+    options_mom_ = set_default_option(options_mom_,'stack_solve_algo',options_.stack_solve_algo);
+    options_mom_.simul = options_.simul;
+    options_mom_.lmmcp = options_.lmmcp;
+    options_mom_.linear_approximation = options_.linear_approximation;
+    options_mom_.endogenous_terminal_period = options_.endogenous_terminal_period;
+    options_mom_.verbosity = options_.verbosity;
+    options_mom_.dynatol = options_.dynatol;
+    options_mom_.mom.initialize_with_previous_endo_simul = false;
+end
+
+% numerical algorithms related
+if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+    options_mom_ = set_default_option(options_mom_,'lyapunov_db',false);                    % doubling algorithm (disclyap_fast) to solve Lyapunov equation to compute variance-covariance matrix of state variables
+    options_mom_ = set_default_option(options_mom_,'lyapunov_fp',false);                    % fixed-point algorithm to solve Lyapunov equation to compute variance-covariance matrix of state variables
+    options_mom_ = set_default_option(options_mom_,'lyapunov_srs',false);                   % square-root-solver (dlyapchol) algorithm to solve Lyapunov equation to compute variance-covariance matrix of state variables
+    options_mom_ = set_default_option(options_mom_,'lyapunov_complex_threshold',1e-15);     % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
+    options_mom_ = set_default_option(options_mom_,'lyapunov_fixed_point_tol',1e-10);       % convergence criterion used in the fixed point Lyapunov solver
+    options_mom_ = set_default_option(options_mom_,'lyapunov_doubling_tol',1e-16);          % convergence criterion used in the doubling algorithm
+    options_mom_ = set_default_option(options_mom_,'sylvester_fp',false);                   % determines whether to use fixed point algorihtm to solve Sylvester equation (gensylv_fp), faster for large scale models
+    options_mom_ = set_default_option(options_mom_,'sylvester_fixed_point_tol',1e-12);      % convergence criterion used in the fixed point Sylvester solver
+end
 
 % Bayesian MCMC related
 if doBayesianEstimation
@@ -332,6 +351,7 @@ options_mom_.block    = options_.block;
 options_mom_.bytecode = options_.bytecode;
 
 % related to steady-state computations
+options_mom_.no_homotopy             = options_.no_homotopy;
 options_mom_.homotopy_force_continue = options_.homotopy_force_continue;
 options_mom_.homotopy_mode           = options_.homotopy_mode;
 options_mom_.homotopy_steps          = options_.homotopy_steps;
@@ -350,8 +370,10 @@ options_mom_.DynareRandomStreams.seed = options_.DynareRandomStreams.seed;
 options_mom_.DynareRandomStreams.algo = options_.DynareRandomStreams.algo;
 
 % dataset_ related
-options_mom_.dataset        = options_.dataset;
-options_mom_.initial_period = options_.initial_period;
+if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+    options_mom_.dataset        = options_.dataset;
+    options_mom_.initial_period = options_.initial_period;
+end
 
 % optimization related
 if any(cellfun(@(x) isnumeric(x) && any(x == 2), options_mom_.optimizer_vec)) % simulated annealing (mode_compute=2)
@@ -404,9 +426,11 @@ options_mom_.analytic_derivation_mode = 0; % needed by get_perturbation_params_d
 options_mom_.initialize_estimated_parameters_with_the_prior_mode = 0; % needed by set_prior.m
 options_mom_.figures = options_.figures; % needed by plot_priors.m
 options_mom_.ramsey_policy = false; % needed by evaluate_steady_state
-options_mom_.risky_steadystate = false;  % needed by resol
-options_mom_.jacobian_flag = true; % needed by dynare_solve
-options_mom_.use_mh_covariance_matrix = false; % needed by posterior_sampler, get's overwritten by similar option in options_mom_.posterior_sampler_options
+if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')    
+    options_mom_.risky_steadystate = false;  % needed by resol
+    options_mom_.jacobian_flag = true; % needed by dynare_solve
+end
+options_mom_.use_mh_covariance_matrix = false; % needed by posterior_sampler, get's overwritten by same option in options_mom_.posterior_sampler_options
 
 % -------------------------------------------------------------------------
 % CHECKS ON SETTINGS
@@ -419,12 +443,14 @@ if strcmp(mom_method,'GMM') || strcmp(mom_method,'SMM')
         error('method_of_moments: Recursive estimation is not supported. Please set an integer as ''first_obs''!');
     end
 end
-if options_mom_.order < 1
-    error('method_of_moments: The order of the Taylor approximation cannot be 0!')
-end
-if options_mom_.order > 2
-    fprintf('Dynare will use ''k_order_solver'' as the order>2\n');
-    options_mom_.k_order_solver = true;
+if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+    if options_mom_.order < 1
+        error('method_of_moments: The order of the Taylor approximation cannot be 0!')
+    end
+    if options_mom_.order > 2
+        fprintf('Dynare will use ''k_order_solver'' as the order>2\n');
+        options_mom_.k_order_solver = true;
+    end
 end
 if strcmp(mom_method,'SMM')
     if options_mom_.mom.simulation_multiple < 1
@@ -433,13 +459,15 @@ if strcmp(mom_method,'SMM')
     end
 end
 if strcmp(mom_method,'GMM')
-    % require pruning with GMM at higher order
-    if options_mom_.order > 1 && ~options_mom_.pruning
-        fprintf('GMM at higher order only works with pruning, so we set pruning option to 1.\n');
-        options_mom_.pruning = true;
-    end
-    if options_mom_.order > 3
-        error('method_of_moments: Perturbation orders higher than 3 are not implemented for GMM estimation, try using SMM!');
+    if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+        % require pruning with GMM at higher order
+        if options_mom_.order > 1 && ~options_mom_.pruning
+            fprintf('GMM at higher order only works with pruning, so we set pruning option to 1.\n');
+            options_mom_.pruning = true;
+        end
+        if options_mom_.order > 3
+            error('method_of_moments: Perturbation orders higher than 3 are not implemented for GMM estimation, try using SMM!');
+        end
     end
 end
 if strcmp(mom_method,'IRF_MATCHING') && doBayesianEstimation
@@ -447,8 +475,8 @@ if strcmp(mom_method,'IRF_MATCHING') && doBayesianEstimation
         warning('method_of_moments: You specified mh_tune_jscale, but the maximum number of iterations is smaller than the step size. No update will take place.')
     end
     if options_mom_.load_results_after_load_mh
-        if ~exist([options_mom_.dirname filesep 'method_of_moments' filesep M_.fname '_mom_results.mat'],'file')
-            fprintf('\nYou specified the ''load_results_after_load_mh'' option, but no ''%s_mom_results.mat'' file\n',M_.fname);
+        if ~exist([options_mom_.dirname filesep 'method_of_moments' filesep fname '_mom_results.mat'],'file')
+            fprintf('\nYou specified the ''load_results_after_load_mh'' option, but no ''%s_mom_results.mat'' file\n',fname);
             fprintf('was found in the folder %s%smethod_of_moments.\n',options_mom_.dirname,filesep);
             fprintf('Results will be recomputed and option ''load_results_after_load_mh'' is reset to false.\n');
             options_mom_.load_results_after_load_mh = false;
diff --git a/matlab/+mom/matched_irfs_blocks.m b/matlab/+mom/matched_irfs_blocks.m
index 64f76ff932db2f5f98a85ae99d7770f084487fb2..1a391552a180a7074e7f82258b7b1876148f905f 100644
--- a/matlab/+mom/matched_irfs_blocks.m
+++ b/matlab/+mom/matched_irfs_blocks.m
@@ -1,5 +1,5 @@
-function [data_irfs, weightMat, irfIndex, maxIrfHorizon] = matched_irfs_blocks(matched_irfs, matched_irfs_weight, varobs_id, obs_nbr, exo_nbr, endo_names)
-% function [data_irfs, weightMat, irfIndex, maxIrfHorizon] = matched_irfs_blocks(matched_irfs, matched_irfs_weight, varobs_id, obs_nbr, exo_nbr, endo_names)
+function [data_irfs, weightMat, irfIndex, maxIrfHorizon] = matched_irfs_blocks(matched_irfs, matched_irfs_weight, varobs_id, obs_nbr, exo_nbr, endo_names, simulation_method)
+% function [data_irfs, weightMat, irfIndex, maxIrfHorizon] = matched_irfs_blocks(matched_irfs, matched_irfs_weight, varobs_id, obs_nbr, exo_nbr, endo_names, simulation_method)
 % -------------------------------------------------------------------------
 % Checks and transforms matched_irfs and matched_irfs_weight blocks
 % for further use in the estimation.
@@ -11,6 +11,7 @@ function [data_irfs, weightMat, irfIndex, maxIrfHorizon] = matched_irfs_blocks(m
 % obs_nbr:             [scalar]     number of observable variables
 % exo_nbr:             [scalar]     number of exogenous variables
 % endo_names:          [cell array] list of endogenous variables
+% simulation_method    [string]     deterministic or stochastic simulation
 % -------------------------------------------------------------------------
 % OUTPUT
 % data_irfs:           [matrix]     irfs for VAROBS as declared in matched_irfs block
@@ -38,43 +39,84 @@ function [data_irfs, weightMat, irfIndex, maxIrfHorizon] = matched_irfs_blocks(m
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 % =========================================================================
-
-maxIrfHorizon = max(cellfun(@(x) x(end), matched_irfs(:,1))); % get maximum irf horizon
-% create full matrix where 1st dimension are irf periods, 2nd dimension are variables as declared in VAROBS, 3rd dimension are shocks.
-data_irfs = nan(maxIrfHorizon,obs_nbr,exo_nbr);
-% overwrite nan values if they are declared in matched_irfs block; remaining nan values will be later ignored in the matching
-for jj = 1:size(matched_irfs,1)
-    idVar       = matched_irfs{jj,1}(1);
-    idVarobs    = find(varobs_id==idVar,1);
-    idShock     = matched_irfs{jj,1}(2);
-    idIrfPeriod = matched_irfs{jj,1}(3);
-    irfValue    = matched_irfs{jj,2};
-    if isempty(idVarobs)
-        skipline;
-        error('method_of_moments: You specified an irf matching involving variable %s, but it is not declared as a varobs!',endo_names{idVar})
+if strcmp(simulation_method,'STOCH_SIMUL')
+    maxIrfHorizon = max(cellfun(@(x) x(end), matched_irfs(:,1))); % get maximum irf horizon
+    % create full matrix where 1st dimension are irf periods, 2nd dimension are variables as declared in VAROBS, 3rd dimension are shocks.
+    data_irfs = nan(maxIrfHorizon,obs_nbr,exo_nbr);
+    % overwrite nan values if they are declared in matched_irfs block; remaining nan values will be later ignored in the matching
+    for jj = 1:size(matched_irfs,1)
+        idVar       = matched_irfs{jj,1}(1);
+        idVarobs    = find(varobs_id==idVar,1);
+        idShock     = matched_irfs{jj,1}(2);
+        idIrfPeriod = matched_irfs{jj,1}(3);
+        irfValue    = matched_irfs{jj,2};
+        if isempty(idVarobs)
+            skipline;
+            error('method_of_moments: You specified an irf matching involving variable %s, but it is not declared as a varobs!',endo_names{idVar})
+        end
+        data_irfs(idIrfPeriod,idVarobs,idShock) = irfValue;
+    end
+    % create (full) empirical weighting matrix
+    weightMat = eye(maxIrfHorizon*obs_nbr*exo_nbr); % identity matrix by default: all irfs are equally important
+    for jj = 1:size(matched_irfs_weight,1)
+        idVar1 = matched_irfs_weight{jj,1}(1);  idVarobs1 = find(varobs_id==idVar1,1);  idShock1 = matched_irfs_weight{jj,1}(2);  idIrfPeriod1 = matched_irfs_weight{jj,1}(3);
+        idVar2 = matched_irfs_weight{jj,2}(1);  idVarobs2 = find(varobs_id==idVar2,1);  idShock2 = matched_irfs_weight{jj,2}(2);  idIrfPeriod2 = matched_irfs_weight{jj,2}(3);
+        weightMatValue = matched_irfs_weight{jj,3};
+        if isempty(idVarobs1)
+            skipline;
+            error('method_of_moments: You specified a weight for an irf matching involving variable %s, but it is not a varobs!',endo_names{idVar1})
+        end
+        if isempty(idVarobs2)s
+            skipline;
+            error('method_of_moments: You specified a weight for an irf matching involving variable %s, but it is not a varobs!',endo_names{idVar2})
+        end
+        idweightMat1 = sub2ind(size(data_irfs),idIrfPeriod1,idVarobs1,idShock1);
+        idweightMat2 = sub2ind(size(data_irfs),idIrfPeriod2,idVarobs2,idShock2);
+        weightMat(idweightMat1,idweightMat2) = weightMatValue;
+        weightMat(idweightMat2,idweightMat1) = weightMatValue; % symmetry
     end
-    data_irfs(idIrfPeriod,idVarobs,idShock) = irfValue;
-end
-% create (full) empirical weighting matrix
-weightMat = eye(maxIrfHorizon*obs_nbr*exo_nbr); % identity matrix by default: all irfs are equally important
-for jj = 1:size(matched_irfs_weight,1)
-    idVar1 = matched_irfs_weight{jj,1}(1);  idVarobs1 = find(varobs_id==idVar1,1);  idShock1 = matched_irfs_weight{jj,1}(2);  idIrfPeriod1 = matched_irfs_weight{jj,1}(3);
-    idVar2 = matched_irfs_weight{jj,2}(1);  idVarobs2 = find(varobs_id==idVar2,1);  idShock2 = matched_irfs_weight{jj,2}(2);  idIrfPeriod2 = matched_irfs_weight{jj,2}(3);
-    weightMatValue = matched_irfs_weight{jj,3};
-    if isempty(idVarobs1)
-        skipline;
-        error('method_of_moments: You specified a weight for an irf matching involving variable %s, but it is not a varobs!',endo_names{idVar1})
+    % focus only on specified irfs
+    irfIndex = find(~isnan(data_irfs));
+    data_irfs = data_irfs(irfIndex);
+    weightMat = weightMat(irfIndex,irfIndex);
+
+elseif strcmp(simulation_method,'PERFECT_FORESIGHT')
+    maxIrfHorizon = max(cellfun(@(x) x(end), matched_irfs(:,1))); % get maximum irf horizon
+    % create full matrix where 1st dimension are irf periods, 2nd dimension are variables as declared in VAROBS.
+    data_irfs = nan(maxIrfHorizon,obs_nbr);
+    % overwrite nan values if they are declared in matched_irfs block; remaining nan values will be later ignored in the matching
+    for jj = 1:size(matched_irfs,1)
+        idVar       = matched_irfs{jj,1}(1);
+        idVarobs    = find(varobs_id==idVar,1);        
+        idIrfPeriod = matched_irfs{jj,1}(2);
+        irfValue    = matched_irfs{jj,2};
+        if isempty(idVarobs)
+            skipline;
+            error('method_of_moments: You specified an irf matching involving variable %s, but it is not declared as a varobs!',endo_names{idVar})
+        end
+        data_irfs(idIrfPeriod,idVarobs) = irfValue;
     end
-    if isempty(idVarobs2)s
-        skipline;
-        error('method_of_moments: You specified a weight for an irf matching involving variable %s, but it is not a varobs!',endo_names{idVar2})
+    % create (full) empirical weighting matrix
+    weightMat = eye(maxIrfHorizon*obs_nbr); % identity matrix by default: all irfs are equally important
+    for jj = 1:size(matched_irfs_weight,1)
+        idVar1 = matched_irfs_weight{jj,1}(1);  idVarobs1 = find(varobs_id==idVar1,1);  idIrfPeriod1 = matched_irfs_weight{jj,1}(2);
+        idVar2 = matched_irfs_weight{jj,2}(1);  idVarobs2 = find(varobs_id==idVar2,1);  idIrfPeriod2 = matched_irfs_weight{jj,2}(2);
+        weightMatValue = matched_irfs_weight{jj,3};
+        if isempty(idVarobs1)
+            skipline;
+            error('method_of_moments: You specified a weight for an irf matching involving variable %s, but it is not a varobs!',endo_names{idVar1})
+        end
+        if isempty(idVarobs2)
+            skipline;
+            error('method_of_moments: You specified a weight for an irf matching involving variable %s, but it is not a varobs!',endo_names{idVar2})
+        end
+        idweightMat1 = sub2ind(size(data_irfs),idIrfPeriod1,idVarobs1);
+        idweightMat2 = sub2ind(size(data_irfs),idIrfPeriod2,idVarobs2);
+        weightMat(idweightMat1,idweightMat2) = weightMatValue;
+        weightMat(idweightMat2,idweightMat1) = weightMatValue; % symmetry
     end
-    idweightMat1 = sub2ind(size(data_irfs),idIrfPeriod1,idVarobs1,idShock1);
-    idweightMat2 = sub2ind(size(data_irfs),idIrfPeriod2,idVarobs2,idShock2);
-    weightMat(idweightMat1,idweightMat2) = weightMatValue;
-    weightMat(idweightMat2,idweightMat1) = weightMatValue; % symmetry
-end
-% focus only on specified irfs
-irfIndex = find(~isnan(data_irfs));
-data_irfs = data_irfs(irfIndex);
-weightMat = weightMat(irfIndex,irfIndex);
\ No newline at end of file
+    % focus only on specified irfs
+    irfIndex = find(~isnan(data_irfs));
+    data_irfs = data_irfs(irfIndex);
+    weightMat = weightMat(irfIndex,irfIndex);
+end
\ No newline at end of file
diff --git a/matlab/+mom/objective_function.m b/matlab/+mom/objective_function.m
index cf9864195364edeb3e9bec587948ac0fa4c644f1..570434adfeb7045d26ed778b2a72e114e6d8ba14 100644
--- a/matlab/+mom/objective_function.m
+++ b/matlab/+mom/objective_function.m
@@ -61,7 +61,11 @@ function [fval, info, exit_flag, df, junkHessian, oo_, M_] = objective_function(
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 % =========================================================================
-
+if options_mom_.mom.initialize_with_previous_endo_simul
+    persistent endo_simul0
+else
+    endo_simul0 = [];
+end
 
 %------------------------------------------------------------------------------
 % Initialization of the returned variables and others...
@@ -99,39 +103,40 @@ end
 
 
 %--------------------------------------------------------------------------
-% Call resol to compute steady state and model solution
+% Stochastic solution: call resol to compute steady state and model solution
 %--------------------------------------------------------------------------
-% Compute linear approximation around the deterministic steady state
-[oo_.dr, info, M_.params] = resol(0, M_, options_mom_, oo_.dr ,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
-% Return, with endogenous penalty when possible, if resol issues an error code
-if info(1)
-    if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 ||...
-            info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1) == 26 || ...
-            info(1) == 81 || info(1) == 84 ||  info(1) == 85 ||  info(1) == 86
-        % meaningful second entry of output that can be used
-        fval = Inf;
-        info(4) = info(2);
-        exit_flag = 0;
-        if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
-            fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
-        end
-        return
-    else
-        fval = Inf;
-        info(4) = 0.1;
-        exit_flag = 0;
-        if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
-            fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
+if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+    % Compute linear approximation around the deterministic steady state
+    [dr, info, M_, oo_] = resol(0, M_, options_mom_, oo_);
+    % Return, with endogenous penalty when possible, if resol issues an error code
+    if info(1)
+        if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 ||...
+                info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1) == 26 || ...
+                info(1) == 81 || info(1) == 84 ||  info(1) == 85 ||  info(1) == 86
+            % meaningful second entry of output that can be used
+            fval = Inf;
+            info(4) = info(2);
+            exit_flag = 0;
+            if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
+                fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
+            end
+            return
+        else
+            fval = Inf;
+            info(4) = 0.1;
+            exit_flag = 0;
+            if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
+                fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
+            end
+            return
         end
-        return
     end
 end
 
-
 %--------------------------------------------------------------------------
-% GMM: Set up pruned state-space system and compute model moments
+% Stochastic GMM: set up pruned state-space system and compute model moments
 %--------------------------------------------------------------------------
-if strcmp(options_mom_.mom.mom_method,'GMM')
+if strcmp(options_mom_.mom.mom_method,'GMM') && strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
     if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
         indpmodel = []; % initialize index for model parameters
         if ~isempty(estim_params_.param_vals)
@@ -210,7 +215,7 @@ end
 
 
 %------------------------------------------------------------------------------
-% SMM: Compute Moments of the model solution for Gaussian innovations
+% Stochastic SMM: compute moments of the model solution for Gaussian innovations
 %------------------------------------------------------------------------------
 if strcmp(options_mom_.mom.mom_method,'SMM')
     % create shock series with correct covariance matrix from iid standard normal shocks
@@ -253,9 +258,8 @@ end
 
 
 %------------------------------------------------------------------------------
-% IRF_MATCHING using STOCH_SIMUL: Compute irfs given model solution and Cholesky
-% decomposition on shock covariance matrix; this resembles the core codes in
-% stoch_simul
+% Stochastic IRF matching: Compute irfs given model solution and Cholesky decomposition
+% on shock covariance matrix; this resembles the core routines in stoch_simul
 %------------------------------------------------------------------------------
 if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') && strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
     cs = get_lower_cholesky_covariance(M_.Sigma_e,options_mom_.add_tiny_number_to_cholesky);
@@ -304,6 +308,86 @@ if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') && strcmp(options_mom_.mom
 end
 
 
+%------------------------------------------------------------------------------
+% Perfect foresight IRF matching: Compute irfs from deterministic simulation
+%------------------------------------------------------------------------------
+if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') && strcmp(options_mom_.mom.simulation_method,'PERFECT_FORESIGHT')
+    noprint_old = options_mom_.noprint;    
+    verbosity_old = options_mom_.verbosity;
+    % oo_.mom.exo_simul = oo_.exo_simul; oo_ = rmfield(oo_,'exo_simul');
+    % oo_.mom.endo_simul = oo_.endo_simul; oo_ = rmfield(oo_,'endo_simul');
+    oo_.mom.exo_steady_state = oo_.exo_steady_state;
+    oo_.mom.exo_det_steady_state = oo_.exo_det_steady_state;
+    % oo_.mom.exo_det_simul = oo_.exo_det_simul; oo_ = rmfield(oo_,'exo_det_simul');
+    oo_.mom.initval_series = oo_.initval_series;
+    oo_.mom.steady_state = oo_.steady_state;
+    options_mom_.noprint = true;
+    [M_, oo_.mom, options_mom_,info] = steady(M_, oo_.mom, options_mom_);
+    options_mom_.noprint = noprint_old;
+    % Return, with endogenous penalty when possible, if steady-state issues an error code
+    if info(1)
+        if info(1) == 19 ||...
+                info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1) == 26 || ...
+                info(1) == 81 || info(1) == 84 ||  info(1) == 85 ||  info(1) == 86
+            % meaningful second entry of output that can be used
+            fval = Inf;
+            info(4) = info(2);
+            exit_flag = 0;
+            if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
+                fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
+            end
+            return
+        else
+            fval = Inf;
+            info(4) = 0.1;
+            exit_flag = 0;
+            if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
+                fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
+            end
+            return
+        end
+    end
+    options_mom_.verbosity = false;
+    options_mom_.noprint = true;
+    oo_.mom = perfect_foresight_setup(M_,options_mom_,oo_.mom);
+    if ~isempty(endo_simul0)
+        oo_.mom.endo_simul = endo_simul0;
+    end
+    oo_.mom = perfect_foresight_solver(M_,options_mom_,oo_.mom,true);
+    if ~oo_.mom.deterministic_simulation.status
+        fval = Inf;
+        info(4) = 0.1;
+        exit_flag = 0;
+        if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
+            fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
+        end
+        return
+    else
+        endo_simul0 = oo_.mom.endo_simul;
+    end
+    options_mom_.noprint = noprint_old;
+    options_mom_.verbosity = verbosity_old;
+    % in deviation from steady-state
+    modelIrf = transpose(oo_.mom.endo_simul(:,2:(options_mom_.irf+1)) - oo_.mom.steady_state);
+    % do transformations on model irfs if irf_matching_file is provided
+    if ~isempty(options_mom_.mom.irf_matching_file.name)
+        [modelIrf,check] = feval(str2func(options_mom_.mom.irf_matching_file.name),modelIrf, M_, options_mom_);
+        if check
+            fval = Inf;
+            info(1) = 180;
+            info(4) = 0.1;
+            exit_flag = 0;
+            if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
+                fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
+            end
+            return
+        end
+    end
+    oo_.mom.irf_model_varobs = modelIrf(:,options_mom_.varobs_id); % focus only on observables (this will be used later for plotting)
+    oo_.mom.model_moments = oo_.mom.irf_model_varobs(oo_.mom.irfIndex); % focus only on selected irf periods
+end
+
+
 %--------------------------------------------------------------------------
 % Compute quadratic target function
 %--------------------------------------------------------------------------
diff --git a/matlab/+mom/print_info_on_estimation_settings.m b/matlab/+mom/print_info_on_estimation_settings.m
index a55a10d2f36183ceef0bfe1a56e8ced2a1ae4a4a..ac479d2b503d1cdc03e9959c3e75005fe036f976 100644
--- a/matlab/+mom/print_info_on_estimation_settings.m
+++ b/matlab/+mom/print_info_on_estimation_settings.m
@@ -114,11 +114,13 @@ for i = 1:length(options_mom_.optimizer_vec)
         fprintf(' (using analytical Jacobian)');
     end
 end
-if options_mom_.order > 0
-    fprintf('\n  - stochastic simulations with perturbation order: %d', options_mom_.order)
-end
-if options_mom_.order > 1 && options_mom_.pruning
-    fprintf(' (with pruning)')
+if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+    if options_mom_.order > 0
+        fprintf('\n  - stochastic simulations with perturbation order: %d', options_mom_.order)
+    end
+    if options_mom_.order > 1 && options_mom_.pruning
+        fprintf(' (with pruning)')
+    end
 end
 if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
     if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors
diff --git a/matlab/+mom/run.m b/matlab/+mom/run.m
index cac0054a131cbd9986fbe1a82aaa4550fc22462b..3e50a0c06aee7ccd7bbe71fcd26a06617dac738a 100644
--- a/matlab/+mom/run.m
+++ b/matlab/+mom/run.m
@@ -143,6 +143,8 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_,
 % - deal with measurement errors (once @wmutschl has implemented this in identification toolbox)
 % - display scaled moments
 % - enable to use first moments even if prefilter option is set
+% SMM
+% - add option to include shocks from other distributions
 % IRF_MATCHING/SMM
 % - add option to do simulations with extended path
 % - add option to do simulations with perfect_foresight and perfect_foresight_with_expectation_errors
@@ -159,7 +161,7 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_,
 % - print more info to console
 
 
-fprintf('\n==== Method of Moments Estimation (%s) ====\n\n',options_mom_.mom.mom_method)
+fprintf('\n==== Method of Moments Estimation (%s) ====\n\n',options_mom_.mom.mom_method);
 
 
 % -------------------------------------------------------------------------
@@ -167,22 +169,25 @@ fprintf('\n==== Method of Moments Estimation (%s) ====\n\n',options_mom_.mom.mom
 % -------------------------------------------------------------------------
 if isempty(estim_params_) % structure storing the info about estimated parameters in the estimated_params block
     if ~(isfield(estim_params_,'nvx') && (size(estim_params_.var_exo,1)+size(estim_params_.var_endo,1)+size(estim_params_.corrx,1)+size(estim_params_.corrn,1)+size(estim_params_.param_vals,1))==0)
-        error('method_of_moments: You need to provide an ''estimated_params'' block!')
+        error('method_of_moments: You need to provide an ''estimated_params'' block!');
     else
-        error('method_of_moments: The ''estimated_params'' block must not be empty!')
+        error('method_of_moments: The ''estimated_params'' block must not be empty!');
     end
 end
 if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
     if ~isfield(M_,'matched_moments') || isempty(M_.matched_moments) % structure storing the moments used for GMM and SMM estimation
-        error('method_of_moments: You need to provide a ''matched_moments'' block for ''mom_method=%s''!',options_mom_.mom.mom_method)
+        error('method_of_moments: You need to provide a ''matched_moments'' block for ''mom_method=%s''!',options_mom_.mom.mom_method);
     end
 elseif strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
     if ~isfield(M_,'matched_irfs') || isempty(M_.matched_irfs) % structure storing the irfs used for matching
-        error('method_of_moments: You need to provide a ''matched_irfs'' block for ''mom_method=%s''!',options_mom_.mom.mom_method)
+        error('method_of_moments: You need to provide a ''matched_irfs'' block for ''mom_method=%s''!',options_mom_.mom.mom_method);
     end
 end
+if (~isempty(estim_params_.var_exo) || ~isempty(estim_params_.var_endo) || ~isempty(estim_params_.corrn) || ~isempty(estim_params_.corrx)) && strcmp(options_mom_.mom.simulation_method, 'PERFECT_FORESIGHT')
+    error('method_of_moments: There is no point in specifying stderr or corr parameters for shocks or measurement error(s) with ''perfect_foresight'' simulations. Please remove these from the ''estimated_params'' block!');
+end
 if (~isempty(estim_params_.var_endo) || ~isempty(estim_params_.corrn)) && strcmp(options_mom_.mom.mom_method, 'GMM')
-    error('method_of_moments: GMM estimation does not support measurement error(s) yet. Please specify them as a structural shock!')
+    error('method_of_moments: GMM estimation does not support measurement error(s) yet. Please specify them as a structural shock!');
 end
 doBayesianEstimation = [estim_params_.var_exo(:,5); estim_params_.var_endo(:,5); estim_params_.corrx(:,6); estim_params_.corrn(:,6); estim_params_.param_vals(:,5)];
 if all(doBayesianEstimation~=0)
@@ -190,10 +195,10 @@ if all(doBayesianEstimation~=0)
 elseif all(doBayesianEstimation==0)
     doBayesianEstimation = false;
 else
-    error('method_of_moments: Estimation must be either fully Frequentist or fully Bayesian. Maybe you forgot to specify a prior distribution!')
+    error('method_of_moments: Estimation must be either fully Frequentist or fully Bayesian. Maybe you forgot to specify a prior distribution!');
 end
 if ~isfield(options_,'varobs')
-    error('method_of_moments: VAROBS statement is missing!')
+    error('method_of_moments: VAROBS statement is missing!');
 end
 check_varobs_are_endo_and_declared_once(options_.varobs,M_.endo_names);
 
@@ -256,11 +261,12 @@ if doBayesianEstimation
 end
 doBayesianEstimationMCMC = doBayesianEstimation && ( (options_mom_.mh_replic>0) || options_mom_.load_mh_file );
 invhess=[];
-% decision rule
-oo_.dr = set_state_space(oo_.dr,M_,options_mom_); % get state-space representation
-oo_.mom.obs_var = []; % create index of observed variables in DR order
-for i = 1:options_mom_.obs_nbr
-    oo_.mom.obs_var = [oo_.mom.obs_var; find(strcmp(options_mom_.varobs{i}, M_.endo_names(oo_.dr.order_var)))];
+if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+    oo_.dr = set_state_space(oo_.dr,M_,options_mom_); % get state-space representation
+    oo_.mom.obs_var = []; % create index of observed variables in DR order
+    for i = 1:options_mom_.obs_nbr
+        oo_.mom.obs_var = [oo_.mom.obs_var; find(strcmp(options_mom_.varobs{i}, M_.endo_names(oo_.dr.order_var)))];
+    end
 end
 
 
@@ -282,7 +288,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
     not_observed_variables=setdiff(oo_.dr.inv_order_var([M_.matched_moments{:,1}]),oo_.mom.obs_var);
     if ~isempty(not_observed_variables)
         skipline;
-        error('method_of_moments: You specified moments involving %s, but it is not a varobs!',M_.endo_names{oo_.dr.order_var(not_observed_variables)})
+        error('method_of_moments: You specified moments involving %s, but it is not a varobs!',M_.endo_names{oo_.dr.order_var(not_observed_variables)});
     end
 end
 
@@ -291,12 +297,13 @@ end
 % matched_irfs: checks and transformations
 % -------------------------------------------------------------------------
 if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
-    [oo_.mom.data_moments, oo_.mom.W, oo_.mom.irfIndex, options_mom_.irf] = mom.matched_irfs_blocks(M_.matched_irfs, M_.matched_irfs_weights, options_mom_.varobs_id, options_mom_.obs_nbr, M_.exo_nbr, M_.endo_names);
+    [oo_.mom.data_moments, oo_.mom.W, oo_.mom.irfIndex, options_mom_.irf] = mom.matched_irfs_blocks(M_.matched_irfs, M_.matched_irfs_weights, options_mom_.varobs_id, options_mom_.obs_nbr, M_.exo_nbr, M_.endo_names, options_mom_.mom.simulation_method);
     oo_.mom.Winv = inv(oo_.mom.W);
     oo_.mom.Winv_logdet = 2*sum(log(diag(chol(oo_.mom.Winv)))); % use this robust option to avoid inf/nan
     options_mom_.mom.mom_nbr = length(oo_.mom.irfIndex);
 end
 
+
 % -------------------------------------------------------------------------
 % irf_matching_file: checks and transformations
 % -------------------------------------------------------------------------
@@ -345,7 +352,9 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
     end
 end
 % check for calibrated covariances before updating parameters
-estim_params_ = check_for_calibrated_covariances(xparam0,estim_params_,M_);
+if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+    estim_params_ = check_for_calibrated_covariances(xparam0,estim_params_,M_);
+end
 % checks on parameter calibration and initialization
 xparam_calib = get_all_parameters(estim_params_,M_); % get calibrated parameters
 if ~any(isnan(xparam_calib)) % all estimated parameters are calibrated
@@ -385,34 +394,40 @@ else
     Bounds.lb = lb;
     Bounds.ub = ub;
     if (strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')) && options_mom_.mom.penalized_estimator
-        fprintf('Penalized estimation turned off as you did not declare priors\n')
+        fprintf('Penalized estimation turned off as you did not declare priors\n');
         options_mom_.mom.penalized_estimator = 0;
     else
         if isfield(options_mom_,'mh_replic') && options_mom_.mh_replic > 0
-            fprintf('Setting ''mh_replic=0'' as you did not declare priors.\n')
+            fprintf('Setting ''mh_replic=0'' as you did not declare priors.\n');
             options_mom_.mh_replic = 0;
         end
     end
 end
 % set correct bounds for standard deviations and correlations
-Bounds = mom.set_correct_bounds_for_stderr_corr(estim_params_,Bounds);
+if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+    Bounds = mom.set_correct_bounds_for_stderr_corr(estim_params_,Bounds);
+end
 % test if initial values of the estimated parameters are all between the prior lower and upper bounds
 if options_mom_.use_calibration_initialization
     try
         check_prior_bounds(xparam0,Bounds,M_,estim_params_,options_mom_,bayestopt_);
     catch last_error
-        fprintf('Cannot use parameter values from calibration as they violate the prior bounds.')
+        fprintf('Cannot use parameter values from calibration as they violate the prior bounds.');
         rethrow(last_error);
     end
 else
     check_prior_bounds(xparam0,Bounds,M_,estim_params_,options_mom_,bayestopt_);
 end
 % check for positive definiteness
-estim_params_ = get_matrix_entries_for_psd_check(M_,estim_params_);
+if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+    estim_params_ = get_matrix_entries_for_psd_check(M_,estim_params_);
+end
 % set sigma_e_is_diagonal flag (needed if the shocks block is not declared in the mod file)
-M_.sigma_e_is_diagonal = true;
-if estim_params_.ncx || any(nnz(tril(M_.Correlation_matrix,-1))) || isfield(estim_params_,'calibrated_covariances')
-    M_.sigma_e_is_diagonal = false;
+if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+    M_.sigma_e_is_diagonal = true;
+    if estim_params_.ncx || any(nnz(tril(M_.Correlation_matrix,-1))) || isfield(estim_params_,'calibrated_covariances')
+        M_.sigma_e_is_diagonal = false;
+    end
 end
 % storing prior parameters in results structure
 if doBayesianEstimation || ( (strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')) && options_mom_.mom.penalized_estimator)
@@ -444,13 +459,13 @@ if doBayesianEstimationMCMC
     end
 end
 % warning if prior allows that stderr parameters are negative or corr parameters are outside the unit circle
-if doBayesianEstimation    
+if doBayesianEstimation && strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
     check_prior_stderr_corr(estim_params_,bayestopt_);
     % check value of prior density
     [~,~,~,info] = priordens(xparam0,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
     if any(info)
-        fprintf('The prior density evaluated at the initial values is Inf for the following parameters: %s\n',bayestopt_.name{info,1})
-        error('The initial value of the prior is -Inf!')
+        fprintf('The prior density evaluated at the initial values is Inf for the following parameters: %s\n',bayestopt_.name{info,1});
+        error('The initial value of the prior is -Inf!');
     end
 end
 
@@ -463,7 +478,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
     % check if datafile has same name as mod file
     [~,name,~] = fileparts(options_mom_.datafile);
     if strcmp(name,M_.fname)
-        error('method_of_moments: ''datafile'' and mod file are not allowed to have the same name; change the name of the ''datafile''!')
+        error('method_of_moments: ''datafile'' and mod file are not allowed to have the same name; change the name of the ''datafile''!');
     end
     dataset_ = makedataset(options_mom_);
     % set options for old interface from the ones for new interface
@@ -479,15 +494,15 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
     % get data moments for the method of moments
     [oo_.mom.data_moments, oo_.mom.m_data] = mom.get_data_moments(dataset_.data, oo_, M_.matched_moments, options_mom_);
     if ~isreal(dataset_.data)
-        error('method_of_moments: The data moments contain complex values!')
+        error('method_of_moments: The data moments contain complex values!');
     end
 end
 
 
 % -------------------------------------------------------------------------
-% SMM: Get shock series and set variance correction factor
+% SMM: get shock series and set variance correction factor
 % -------------------------------------------------------------------------
-if strcmp(options_mom_.mom.mom_method,'SMM')
+if strcmp(options_mom_.mom.mom_method,'SMM') && strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
     options_mom_.mom.long = round(options_mom_.mom.simulation_multiple*options_mom_.nobs);
     options_mom_.mom.variance_correction_factor = (1+1/options_mom_.mom.simulation_multiple);
     % draw shocks for SMM
@@ -527,7 +542,7 @@ else
 end
 [oo_.steady_state, info, steady_state_changes_parameters] = check_steady_state_changes_parameters(M_, estim_params_, oo_, options_mom_, steadystate_check_flag_vec);
 if info(1)
-    fprintf('\nThe steady state at the initial parameters cannot be computed.\n')
+    fprintf('\nThe steady state at the initial parameters cannot be computed.\n');
     print_info(info, 0, options_mom_);
 end
 if steady_state_changes_parameters && strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors
@@ -539,6 +554,19 @@ end
 test_for_deep_parameters_calibration(M_);
 
 
+% -------------------------------------------------------------------------
+% Perfect foresight: get declared blocks
+% -------------------------------------------------------------------------
+if strcmp(options_mom_.mom.simulation_method,'PERFECT_FORESIGHT')
+    % oo_.mom.exo_simul = oo_.exo_simul; oo_ = rmfield(oo_,'exo_simul');
+    % oo_.mom.endo_simul = oo_.endo_simul; oo_ = rmfield(oo_,'endo_simul');
+    % oo_.mom.exo_steady_state = oo_.exo_steady_state;
+    % oo_.mom.exo_det_steady_state = oo_.exo_det_steady_state;
+    % oo_.mom.exo_det_simul = oo_.exo_det_simul; oo_ = rmfield(oo_,'exo_det_simul');
+    % oo_.mom.initval_series = oo_.initval_series; oo_ = rmfield(oo_,'initval_series');
+    % oo_.mom.steady_state = oo_.steady_state;
+end
+
 % -------------------------------------------------------------------------
 % checks for objective function at initial parameters
 % -------------------------------------------------------------------------
@@ -553,20 +581,26 @@ try
     elapsed_time = toc(tic_id);
     if isnan(fval)
         if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
-            error('method_of_moments: The initial value of the objective function with identity weighting matrix is NaN!')
+            error('method_of_moments: The initial value of the objective function with identity weighting matrix is NaN!');
         else
-            error('method_of_moments: The initial value of the objective function is NaN!')
+            error('method_of_moments: The initial value of the objective function is NaN!');
         end
     elseif imag(fval)
         if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
-            error('method_of_moments: The initial value of the objective function with identity weighting matrix is complex!')
+            error('method_of_moments: The initial value of the objective function with identity weighting matrix is complex!');
+        else
+            error('method_of_moments: The initial value of the objective function is complex!');
+        end
+    elseif isinf(fval)
+        if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
+            error('method_of_moments: The initial value of the objective function with identity weighting matrix is infinity!');
         else
-            error('method_of_moments: The initial value of the objective function is complex!')
+            error('method_of_moments: The initial value of the objective function is infinity!');
         end
     end
     if info(1) > 0
-        disp('method_of_moments: Error in computing the objective function for initial parameter values')
-        print_info(info, options_mom_.noprint, options_mom_)
+        disp('method_of_moments: Error in computing the objective function for initial parameter values');
+        print_info(info, options_mom_.noprint, options_mom_);
     end
     fprintf('Initial value of the moment objective function');
     if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
@@ -577,10 +611,10 @@ try
 catch last_error % if check fails, provide info on using calibration if present
     if estim_params_.full_calibration_detected %calibrated model present and no explicit starting values
         skipline(1);
-        fprintf('There was an error in computing the moments for initial parameter values.\n')
-        fprintf('If this is not a problem with the setting of options (check the error message below),\n')
-        fprintf('you should try using the calibrated version of the model as starting values. To do\n')
-        fprintf('this, add an empty estimated_params_init-block with use_calibration option immediately before the estimation\n')
+        fprintf('There was an error in computing the moments for initial parameter values.\n');
+        fprintf('If this is not a problem with the setting of options (check the error message below),\n');
+        fprintf('you should try using the calibrated version of the model as starting values. To do\n');
+        fprintf('this, add an empty estimated_params_init-block with use_calibration option immediately before the estimation\n');
         fprintf('command (and after the estimated_params-block so that it does not get overwritten):\n');
         skipline(2);
     end
@@ -776,13 +810,13 @@ if doBayesianEstimationMCMC
             % THIS IS PROBABLY NOT USEFUL HERE AND CAN BE REMOVED (PREPROCESSOR: REMOVE bayesian_irf, moments_varendo)
             %if options_mom_.bayesian_irf
             %    if error_flag
-            %        error('method_of_moments: Cannot compute the posterior IRFs!')
+            %        error('method_of_moments: Cannot compute the posterior IRFs!');
             %    end
             %    PosteriorIRF('posterior','method_of_moments::mcmc');
             %end
             % if options_mom_.moments_varendo
             %     if error_flag
-            %         error('method_of_moments: Cannot compute the posterior moments for the endogenous variables!')
+            %         error('method_of_moments: Cannot compute the posterior moments for the endogenous variables!');
             %     end
             %     if options_mom_.load_mh_file && options_mom_.mh_replic==0 %user wants to recompute results
             %        [MetropolisFolder, info] = CheckPath('metropolis',options_mom_.dirname);
@@ -806,7 +840,7 @@ if doBayesianEstimationMCMC
             %     oo_ = compute_moments_varendo('posterior',options_,M_,oo_,var_list_);
             % end            
         else
-            fprintf('''sub_draws'' was set to 0. Skipping posterior computations.')
+            fprintf('''sub_draws'' was set to 0. Skipping posterior computations.');
         end
         xparam1 = get_posterior_parameters('mean',M_,estim_params_,oo_.mom,options_);
     end
@@ -825,7 +859,7 @@ if doBayesianEstimationMCMC
     % end
     % disp(' ');
     % disp(['Correlations of Parameters (at Posterior Mode) > ',num2str(ExtremeCorrBound)]);
-    % disp(ExtremeCorrParams)
+    % disp(ExtremeCorrParams);
 end
 
 
@@ -839,15 +873,19 @@ if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_meth
     oo_ = mom.Jtest(xparam1, objective_function, Woptflag, oo_, options_mom_, bayestopt_, Bounds, estim_params_, M_, dataset_.nobs);    
 elseif strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
     if ~options_mom_.nograph
-        mom.graph_comparison_irfs(M_.matched_irfs,oo_.mom.irf_model_varobs,options_mom_.varobs_id,options_mom_.irf,options_mom_.relative_irf,M_.endo_names,M_.exo_names,M_.exo_names_tex,M_.dname,M_.fname,options_mom_.graph_format,options_mom_.TeX,options_mom_.nodisplay,options_mom_.figures.textwidth)
+        if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+            mom.graph_comparison_irfs(M_.matched_irfs,oo_.mom.irf_model_varobs,options_mom_.varobs_id,options_mom_.irf,options_mom_.relative_irf,M_.endo_names,M_.exo_names,M_.exo_names_tex,M_.dname,M_.fname,options_mom_.graph_format,options_mom_.TeX,options_mom_.nodisplay,options_mom_.figures.textwidth)
+        end
     end
 end
 % display comparison of model moments/irfs and data moments/irfs
-mom.display_comparison_moments_irfs(M_, options_mom_, oo_.mom.data_moments, oo_.mom.model_moments);
+if strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
+    mom.display_comparison_moments_irfs(M_, options_mom_, oo_.mom.data_moments, oo_.mom.model_moments);
+end
 % save results to _mom_results.mat
 save([M_.dname filesep 'method_of_moments' filesep 'cet_rwmh_mom_results.mat'], 'oo_', 'options_mom_', 'M_', 'estim_params_', 'bayestopt_'); 
 
-fprintf('\n==== Method of Moments Estimation (%s) Completed ====\n\n',options_mom_.mom.mom_method)
+fprintf('\n==== Method of Moments Estimation (%s) Completed ====\n\n',options_mom_.mom.mom_method);
 
 
 % -------------------------------------------------------------------------
diff --git a/matlab/perfect-foresight-models/perfect_foresight_setup.m b/matlab/perfect-foresight-models/perfect_foresight_setup.m
index f187b9dc2e38ddd6e7a7449b581c592d7d20c527..89c8c2a293e2dff25b25e349ffd9de22e1d3c8e9 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_setup.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_setup.m
@@ -1,11 +1,18 @@
-function perfect_foresight_setup()
+function oo_ = perfect_foresight_setup(M_,options_,oo_)
 % Prepares a deterministic simulation, by filling oo_.exo_simul and oo_.endo_simul
 %
 % INPUTS
-%   None
+%  M_          [structure] model structure
+%  options_    [structure] options structure
+%  oo_         [structure] output structure
 %
 % OUTPUTS
-%   none
+% oo_.endo_simul            [matrix]  guess values for endogenous variables
+% oo_.exo_det_simul         [matrix]  scenario values for deterministic exogenous variables
+% oo_.exo_det_steady_state  [matrix]  scenario values for deterministic exogenous variables in steady-state (if previously empty)
+% oo_.exo_simul             [matrix]  scenario values for exogenous variables
+% oo_.exo_steady_state      [matrix]  scenario values for exogenous variables in steady-state (if previously empty)
+% oo_.steady_state          [matrix]  guess values for endogenous variables in steady-state (if previously empty)
 %
 % ALGORITHM
 %
@@ -29,7 +36,9 @@ function perfect_foresight_setup()
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-global M_ options_ oo_
+if nargin == 0
+    global M_ options_ oo_
+end
 
 test_for_deep_parameters_calibration(M_);
 
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m
index 97f53848fd5e52a42a340f7edff1e4b0175c6c87..bbc77bddb3f340654a9232644e06ae07470f607a 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m
@@ -1,4 +1,4 @@
-function perfect_foresight_solver()
+function oo_ = perfect_foresight_solver(M_, options_, oo_, isEstimation)
 % Computes deterministic simulations
 %
 % INPUTS
@@ -28,8 +28,10 @@ function perfect_foresight_solver()
 %
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
-
-global M_ options_ oo_ ys0_ ex0_
+global ys0_ ex0_
+if nargin == 0
+    global M_ options_ oo_
+end
 
 check_input_arguments(options_, M_, oo_);
 
@@ -302,7 +304,7 @@ if ~isreal(oo_.endo_simul(:)) % cannot happen with bytecode or the perfect_fores
         oo_.endo_simul = real_simul;
         maxerror = real_maxerror;
     else
-        current_share = 0;
+        current_share = 0;        
         disp('Simulation terminated with imaginary parts in the residuals or endogenous variables.')
     end
 end
@@ -376,27 +378,29 @@ if did_homotopy
     warning(warning_old_state);
 end
 
-dyn2vec(M_, oo_, options_);
-
-if isfield(oo_, 'initval_series') && ~isempty(oo_.initval_series)
-    initial_period = oo_.initval_series.dates(1)+(M_.orig_maximum_lag-1);
-elseif ~isdates(options_.initial_period) && isnan(options_.initial_period)
-    initial_period = dates(1,1);
-else
-    initial_period = options_.initial_period;
-end
-
-ts = dseries([transpose(oo_.endo_simul(1:M_.orig_endo_nbr,:)), oo_.exo_simul], initial_period, [M_.endo_names(1:M_.orig_endo_nbr); M_.exo_names]);
-
-if isfield(oo_, 'initval_series') && ~isempty(oo_.initval_series)
-    names = ts.name;
-    ts = merge(oo_.initval_series{names{:}}, ts);
-end
-
-assignin('base', 'Simulated_time_series', ts);
-
-if success
-    oo_.gui.ran_perfect_foresight = true;
+if ~isEstimation
+    dyn2vec(M_, oo_, options_);
+    
+    if isfield(oo_, 'initval_series') && ~isempty(oo_.initval_series)
+        initial_period = oo_.initval_series.dates(1)+(M_.orig_maximum_lag-1);
+    elseif ~isdates(options_.initial_period) && isnan(options_.initial_period)
+        initial_period = dates(1,1);
+    else
+        initial_period = options_.initial_period;
+    end
+    
+    ts = dseries([transpose(oo_.endo_simul(1:M_.orig_endo_nbr,:)), oo_.exo_simul], initial_period, [M_.endo_names(1:M_.orig_endo_nbr); M_.exo_names]);
+    
+    if isfield(oo_, 'initval_series') && ~isempty(oo_.initval_series)
+        names = ts.name;
+        ts = merge(oo_.initval_series{names{:}}, ts);
+    end
+    
+    assignin('base', 'Simulated_time_series', ts);
+    
+    if success
+        oo_.gui.ran_perfect_foresight = true;
+    end
 end
 
 end
diff --git a/matlab/steady.m b/matlab/steady.m
index 255f4d7db944eb18fbc5d39369af1c58b3e26bd2..6b2f2480fb5e258dfeca8f84c004ee709679ce87 100644
--- a/matlab/steady.m
+++ b/matlab/steady.m
@@ -1,4 +1,4 @@
-function steady()
+function [M_, oo_, options_,info] = steady(M_, oo_, options_)
 % function steady()
 % computes and prints the steady state calculations
 %
@@ -28,7 +28,9 @@ function steady()
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-global M_ oo_ options_
+if nargin == 0
+    global M_ oo_ options_
+end
 
 test_for_deep_parameters_calibration(M_);