diff --git a/matlab/PosteriorIRF_core1.m b/matlab/PosteriorIRF_core1.m
index 82c5cfab9ada0d5bf8dd380ca81822779b3376ba..403c0b557aa2e98472d58cb25c77a79f9d0b5814 100644
--- a/matlab/PosteriorIRF_core1.m
+++ b/matlab/PosteriorIRF_core1.m
@@ -176,9 +176,9 @@ while fpar<B
     for i=irf_shocks_indx
         if SS(i,i) > 1e-13
             if options_.order>1 && options_.relative_irf % normalize shock to 0.01 before IRF generation for GIRFs; multiply with 100 later
-                y=irf(dr,SS(M_.exo_names_orig_ord,i)./SS(i,i)/100, options_.irf, options_.drop,options_.replic,options_.order);
+                y=irf(M_,options_,dr,SS(M_.exo_names_orig_ord,i)./SS(i,i)/100, options_.irf, options_.drop,options_.replic,options_.order);
             else
-                y=irf(dr,SS(M_.exo_names_orig_ord,i), options_.irf, options_.drop,options_.replic,options_.order);
+                y=irf(M_,options_,dr,SS(M_.exo_names_orig_ord,i), options_.irf, options_.drop,options_.replic,options_.order);
             end
             if options_.relative_irf && options_.order==1 %multiply with 100 for backward compatibility
                 y = 100*y/SS(i,i);
diff --git a/matlab/discretionary_policy.m b/matlab/discretionary_policy.m
index c565353ce7813310db71200150f2ff00f27d7479..521a6b351bffae3db80f40b062545f0183380ef6 100644
--- a/matlab/discretionary_policy.m
+++ b/matlab/discretionary_policy.m
@@ -1,6 +1,6 @@
-function info = discretionary_policy(var_list)
+function [info, oo_, options_] = discretionary_policy(M_, options_, oo_, var_list)
 
-% Copyright (C) 2007-2018 Dynare Team
+% Copyright (C) 2007-2019 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -17,12 +17,10 @@ function info = discretionary_policy(var_list)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global options_ oo_ M_
-
+origorder = options_.order;
 options_.discretionary_policy = 1;
-oldoptions = options_;
 options_.order = 1;
-info = stoch_simul(var_list);
+[info, oo_] = stoch_simul(M_, options_, oo_, var_list);
 
 if ~options_.noprint
     disp_steady_state(M_,oo_)
@@ -35,4 +33,4 @@ end
 
 oo_.planner_objective_value = evaluate_planner_objective(M_,options_,oo_);
 
-options_ = oldoptions;
+options_.order = origorder;
diff --git a/matlab/disp_moments.m b/matlab/disp_moments.m
index ee62fb3f7e6dada898c15c0cb5e705f5905c8136..6ae4e4d33997a0255ab5aad68a26a5a1205537e9 100644
--- a/matlab/disp_moments.m
+++ b/matlab/disp_moments.m
@@ -11,7 +11,7 @@ function oo_=disp_moments(y,var_list,M_,options_,oo_)
 % OUTPUTS
 %   oo_                 [structure]    Dynare's results structure,
 
-% Copyright (C) 2001-2018 Dynare Team
+% Copyright (C) 2001-2019 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -166,7 +166,7 @@ if ~options_.nodecomposition
             temp_shock_mat=zeros(size(shock_mat));
             temp_shock_mat(:,i_exo_var(shock_iter))=shock_mat(:,i_exo_var(shock_iter));
             temp_shock_mat(:,i_exo_var) = temp_shock_mat(:,i_exo_var)*chol_S;
-            y_sim_one_shock = simult_(y0,oo_.dr,temp_shock_mat,options_.order);
+            y_sim_one_shock = simult_(M_,options_,y0,oo_.dr,temp_shock_mat,options_.order);
             y_sim_one_shock=y_sim_one_shock(ivar,1+options_.drop+1:end)';
             y_sim_one_shock=get_filtered_time_series(y_sim_one_shock,mean(y_sim_one_shock),options_);
             oo_.variance_decomposition(:,i_exo_var(shock_iter))=var(y_sim_one_shock)./s2*100;            
diff --git a/matlab/dyn2vec.m b/matlab/dyn2vec.m
index aabb976b1ce390f3d818b6c01371e9d492ac8a03..46a3b5adfcbc8f907f9378bbb7017890ab959a3f 100644
--- a/matlab/dyn2vec.m
+++ b/matlab/dyn2vec.m
@@ -1,5 +1,5 @@
-function [z,zss]=dyn2vec(s1,s2)
-% function [z,zss]=dyn2vec(s1,s2)
+function [z,zss]=dyn2vec(M_, oo_, options_, s1, s2)
+% function [z,zss]=dyn2vec(M_, oo_, options_, s1, s2)
 % Takes Dynare variables from oo_.endo_simul and copies them into matlab global vectors
 %
 % INPUTS
@@ -14,7 +14,7 @@ function [z,zss]=dyn2vec(s1,s2)
 %   none
 %
 
-% Copyright (C) 2001-2018 Dynare Team
+% Copyright (C) 2001-2019 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -31,7 +31,9 @@ function [z,zss]=dyn2vec(s1,s2)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global M_ oo_ options_
+if ~(nargin >= 3)
+    error('DYNARE dyn2vec error: function takes at least 3 arguments');
+end
 
 if options_.smpl == 0
     k = [1:size(oo_.endo_simul,2)];
@@ -39,7 +41,7 @@ else
     k = [M_.maximum_lag+options_.smpl(1):M_.maximum_lag+options_.smpl(2)];
 end
 
-if nargin == 0
+if nargin == 3
     if nargout > 0
         t = ['DYNARE dyn2vec error: the function doesn''t return values when' ...
              ' used without input argument'];
@@ -70,9 +72,9 @@ else
 end
 
 if nargout == 0
-    if nargin == 1
+    if nargin == 4
         assignin('base',s1,z);
-    elseif nargin == 2
+    elseif nargin == 5
         assignin('base',s2,z);
     end
 else
diff --git a/matlab/ep/extended_path_core.m b/matlab/ep/extended_path_core.m
index 5a99b6f21d897bc90a58c68ee4d6d4392b6b22be..1768c0bbd60731f0813270540c766c37c6497d4b 100644
--- a/matlab/ep/extended_path_core.m
+++ b/matlab/ep/extended_path_core.m
@@ -4,7 +4,7 @@ function [y, info_convergence, endogenousvariablespaths] = extended_path_core(pe
                                                   debug,bytecode_flag,order,M,pfm,algo,solve_algo,stack_solve_algo,...
                                                   olmmcp,options,oo,initialguess)
 
-% Copyright (C) 2016-2017 Dynare Team
+% Copyright (C) 2016-2019 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -24,7 +24,7 @@ function [y, info_convergence, endogenousvariablespaths] = extended_path_core(pe
 ep = options.ep;
 
 if init% Compute first order solution (Perturbation)...
-    endo_simul = simult_(initial_conditions,oo.dr,exo_simul(2:end,:),1);
+    endo_simul = simult_(M,options,initial_conditions,oo.dr,exo_simul(2:end,:),1);
 else
     if nargin==20 && ~isempty(initialguess)
         % Note that the first column of initialguess should be equal to initial_conditions.
diff --git a/matlab/forcst.m b/matlab/forcst.m
index 67bc0f7e2192b5a0bf295dc46889935aa4c51edf..4aa28230b6420cfb10a037f6c19d633ffd23fbab 100644
--- a/matlab/forcst.m
+++ b/matlab/forcst.m
@@ -22,7 +22,7 @@ function [yf,int_width,int_width_ME]=forcst(dr,y0,horizon,var_list,M_,oo_,option
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2003-2018 Dynare Team
+% Copyright (C) 2003-2019 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -40,7 +40,7 @@ function [yf,int_width,int_width_ME]=forcst(dr,y0,horizon,var_list,M_,oo_,option
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 oo_=make_ex_(M_,options_,oo_);
-yf = simult_(y0,dr,zeros(horizon,M_.exo_nbr),1);
+yf = simult_(M_,options_,y0,dr,zeros(horizon,M_.exo_nbr),1);
 nstatic = M_.nstatic;
 nspred = M_.nspred;
 nc = size(dr.ghx,2);
diff --git a/matlab/gsa/stab_map_.m b/matlab/gsa/stab_map_.m
index 15c6f22a5ca21b439979bfc187ffd5be805cf1ea..6758c41fe21c14d56721e15fc712313a535c35c4 100644
--- a/matlab/gsa/stab_map_.m
+++ b/matlab/gsa/stab_map_.m
@@ -495,7 +495,7 @@ else
         options_.nomoments=1;
         options_.irf=0;
         options_.noprint=1;
-        stoch_simul([]);
+        [~, oo_, options_] = stoch_simul(M_, options_, oo_, []);
         %T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),length(istable));
         ntrans=length(istable);
         yys=NaN(length(ys_),ntrans);
diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m
index 9939a1590e26712588febce74501d05aafb5ca1a..9f72aaafb8631027fd6ef40b2e8b71ca66d387ed 100644
--- a/matlab/identification_analysis.m
+++ b/matlab/identification_analysis.m
@@ -303,7 +303,7 @@ if info(1) == 0 %no errors in solution
                 end
                 analytic_derivation              = options_.analytic_derivation;
                 options_.analytic_derivation     = -2; %this sets asy_Hess=1 in dsge_likelihood.m                
-                info = stoch_simul(options_.varobs);
+                [info, oo_, options_] = stoch_simul(M_, options_, oo_, options_.varobs);
                 dataset_ = dseries(oo_.endo_simul(options_.varobs_id,100+1:end)',dates('1Q1'), options_.varobs); %get information on moments
                 derivatives_info.no_DLIK = 1;
                 bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding
diff --git a/matlab/irf.m b/matlab/irf.m
index 2355392d76c32bfe1a90491208104e856f9c4cbb..d1c7014e87c25e46aa75dbad6279759a2f009269 100644
--- a/matlab/irf.m
+++ b/matlab/irf.m
@@ -1,15 +1,16 @@
-function y = irf(dr, e1, long, drop, replic, iorder)
-
-% function y = irf(dr, e1, long, drop, replic, iorder)
+function y = irf(M_, options_, dr, e1, long, drop, replic, iorder)
+% function y = irf(M_, options_, dr, e1, long, drop, replic, iorder)
 % Computes impulse response functions
 %
 % INPUTS
-%    dr:     structure of decisions rules for stochastic simulations
-%    e1:     exogenous variables value in time 1 after one shock
-%    long:   number of periods of simulation
-%    drop:   truncation (in order 2)
-%    replic: number of replications (in order 2)
-%    iorder: first or second order approximation
+%    M_:        structure representing the model
+%    options_:  structure representing options for commands
+%    dr:        structure of decisions rules for stochastic simulations
+%    e1:        exogenous variables value in time 1 after one shock
+%    long:      number of periods of simulation
+%    drop:      truncation (in order 2)
+%    replic:    number of replications (in order 2)
+%    iorder:    first or second order approximation
 %
 % OUTPUTS
 %    y:      impulse response matrix
@@ -17,7 +18,7 @@ function y = irf(dr, e1, long, drop, replic, iorder)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2003-2017 Dynare Team
+% Copyright (C) 2003-2019 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -34,9 +35,6 @@ function y = irf(dr, e1, long, drop, replic, iorder)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global M_ oo_ options_
-
-
 if M_.maximum_lag >= 1
     temps = repmat(dr.ys,1,M_.maximum_lag);
 else
@@ -53,21 +51,20 @@ if local_order == 1
     y1 = repmat(dr.ys,1,long);
     ex2 = zeros(long,M_.exo_nbr);
     ex2(1,:) = e1';
-    y2 = simult_(temps,dr,ex2,local_order);
+    y2 = simult_(M_,options_,temps,dr,ex2,local_order);
     y = y2(:,M_.maximum_lag+1:end)-y1;
 else
     % eliminate shocks with 0 variance
     i_exo_var = setdiff([1:M_.exo_nbr],find(diag(M_.Sigma_e) == 0 ));
     nxs = length(i_exo_var);
     ex1 = zeros(long+drop,M_.exo_nbr);
-    ex2 = ex1;
     chol_S = chol(M_.Sigma_e(i_exo_var,i_exo_var));
     for j = 1: replic
         ex1(:,i_exo_var) = randn(long+drop,nxs)*chol_S;
         ex2 = ex1;
         ex2(drop+1,:) = ex2(drop+1,:)+e1';
-        y1 = simult_(temps,dr,ex1,local_order);
-        y2 = simult_(temps,dr,ex2,local_order);
+        y1 = simult_(M_,options_,temps,dr,ex1,local_order);
+        y2 = simult_(M_,options_,temps,dr,ex2,local_order);
         y = y+(y2(:,M_.maximum_lag+drop+1:end)-y1(:,M_.maximum_lag+drop+1:end));
     end
     y=y/replic;
diff --git a/matlab/osr.m b/matlab/osr.m
index d6d9786c59948b66eda1bf03bdec04082f7e5634..25db19bdd4e52308f30535908607376d589dc131 100644
--- a/matlab/osr.m
+++ b/matlab/osr.m
@@ -62,4 +62,4 @@ if ~options_.noprint
 end
 osr_res = osr1(i_params,i_var,W);
 
-stoch_simul(var_list);
\ No newline at end of file
+[~, oo_, options_] = stoch_simul(M_, options_, oo_, var_list);
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m
index 4c55c28232ccfcc97f0e46fe6ff59ce1db6af9b4..8a90e6000581382d4079e6ffd04a7fba4af500fc 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m
@@ -208,7 +208,7 @@ end
 
 skipline()
 
-dyn2vec;
+dyn2vec(M_, oo_, options_);
 
 if ~isdates(options_.initial_period) && isnan(options_.initial_period)
     initial_period = dates(1,1);
diff --git a/matlab/ramsey_policy.m b/matlab/ramsey_policy.m
index 145aa56e6b61c53e81b67c91c343e7dc74d70db1..0ce386c8a3c7d0ae92d001187f3f0fadb5c932d9 100644
--- a/matlab/ramsey_policy.m
+++ b/matlab/ramsey_policy.m
@@ -1,6 +1,6 @@
 function info = ramsey_policy(var_list)
 
-% Copyright (C) 2007-2018 Dynare Team
+% Copyright (C) 2007-2019 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -39,7 +39,7 @@ else
     end
 end
 
-info = stoch_simul(var_list);
+[info, oo_, options_] = stoch_simul(M_, options_, oo_, var_list);
 
 oo_.steady_state = oo_.dr.ys;
 
diff --git a/matlab/simult.m b/matlab/simult.m
index 5c7dd9b7c9fbabffa7e0a4de2bde43b6f6cf1fb8..f96827f75b2f727448934e7c4d79f4fc1a57ae4d 100644
--- a/matlab/simult.m
+++ b/matlab/simult.m
@@ -45,7 +45,7 @@ function [y_out,DynareResults] =simult(y0, dr,DynareModel,DynareOptions,DynareRe
 %! @end deftypefn
 %@eod:
 
-% Copyright (C) 2001-2012 Dynare Team
+% Copyright (C) 2001-2019 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -82,7 +82,7 @@ for i=1:replic
         % independently of the length of the simulation
         DynareResults.exo_simul(:,i_exo_var) = randn(nxs,DynareOptions.periods)'*chol_S;
     end
-    y_ = simult_(y0,dr,DynareResults.exo_simul,order);
+    y_ = simult_(DynareModel,DynareOptions,y0,dr,DynareResults.exo_simul,order);
     % elimninating initial value
     y_ = y_(:,2:end);
     if replic > 1
diff --git a/matlab/simult_.m b/matlab/simult_.m
index e4fa39618d78d5d464e60893062f8811a5264614..e3ec1ec31724d2e9d4bd9816a2f6b38c2a1e9af2 100644
--- a/matlab/simult_.m
+++ b/matlab/simult_.m
@@ -1,8 +1,10 @@
-function y_=simult_(y0,dr,ex_,iorder)
+function y_=simult_(M_,options_,y0,dr,ex_,iorder)
 % Simulates the model using a perturbation approach, given the path for the exogenous variables and the
 % decision rules.
 %
 % INPUTS
+%    M_       [struct]   model
+%    options_ [struct]   options
 %    y0       [double]   n*1 vector, initial value (n is the number of declared endogenous variables plus the number
 %                        of auxilliary variables for lags and leads); must be in declaration order, i.e. as in M_.endo_names
 %    dr       [struct]   matlab's structure where the reduced form solution of the model is stored.
@@ -15,7 +17,7 @@ function y_=simult_(y0,dr,ex_,iorder)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2001-2017 Dynare Team
+% Copyright (C) 2001-2019 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -32,8 +34,6 @@ function y_=simult_(y0,dr,ex_,iorder)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global M_ options_
-
 iter = size(ex_,1);
 endo_nbr = M_.endo_nbr;
 exo_nbr = M_.exo_nbr;
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index 819638f1a64393028143f1f34d0e819f1900a21f..9efb10530b3121f737973a843b73dd4c3ccb830a 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -1,4 +1,4 @@
-function info = stoch_simul(var_list)
+function [info, oo_, options_] = stoch_simul(M_, options_, oo_, var_list)
 
 % Copyright (C) 2001-2019 Dynare Team
 %
@@ -17,8 +17,6 @@ function info = stoch_simul(var_list)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global M_ options_ oo_
-
 % Test if the order of approximation is nonzero (the preprocessor tests if order is non negative).
 if isequal(options_.order,0)
     error('stoch_simul:: The order of the Taylor approximation cannot be 0!')
@@ -90,7 +88,7 @@ else
         oo_.steady_state=exp(oo_.steady_state);
         options_.logged_steady_state=0;
     end
-    [oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
+    [~,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
 end
 
 if options_.loglinear && isfield(oo_.dr,'ys') && options_.logged_steady_state==0 %log steady state for correct display of decision rule
@@ -182,7 +180,7 @@ if options_.periods > 0 && ~PI_PCL_solver
     [ys, oo_] = simult(y0,oo_.dr,M_,options_,oo_);
     oo_.endo_simul = ys;
     if ~options_.minimal_workspace
-        dyn2vec;
+        dyn2vec(M_, oo_, options_);
     end
 end
 
@@ -222,10 +220,10 @@ if options_.irf
                 y=PCL_Part_info_irf (0, PCL_varobs, i_var, M_, oo_.dr, options_.irf, i);
             else
                 if options_.order>1 && options_.relative_irf % normalize shock to 0.01 before IRF generation for GIRFs; multiply with 100 later
-                    y=irf(oo_.dr,cs(M_.exo_names_orig_ord,i)./cs(i,i)/100, options_.irf, options_.drop, ...
+                    y=irf(M_, options_, oo_.dr,cs(M_.exo_names_orig_ord,i)./cs(i,i)/100, options_.irf, options_.drop, ...
                           options_.replic, options_.order);
                 else %for linear model, rescaling is done later
-                    y=irf(oo_.dr,cs(M_.exo_names_orig_ord,i), options_.irf, options_.drop, ...
+                    y=irf(M_, options_, oo_.dr,cs(M_.exo_names_orig_ord,i), options_.irf, options_.drop, ...
                           options_.replic, options_.order);
                 end
             end
diff --git a/preprocessor b/preprocessor
index bad0c3cf271fd4a9ed7729096f32bd078e59f486..04b6a68aef1dea23f991e788af911781d41bb4ab 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit bad0c3cf271fd4a9ed7729096f32bd078e59f486
+Subproject commit 04b6a68aef1dea23f991e788af911781d41bb4ab
diff --git a/tests/conditional_forecasts/3/fs2000_conditional_forecast_initval.mod b/tests/conditional_forecasts/3/fs2000_conditional_forecast_initval.mod
index f37287e2a0c06a4406be8586425f9e1eb93bf48e..95481d2dc834b043564e2c6bdc0b6174672658de 100644
--- a/tests/conditional_forecasts/3/fs2000_conditional_forecast_initval.mod
+++ b/tests/conditional_forecasts/3/fs2000_conditional_forecast_initval.mod
@@ -97,7 +97,7 @@ shock_matrix = zeros(options_cond_fcst_.periods ,M_.exo_nbr); %create shock matr
 shock_matrix(1:5,strmatch('e_a',M_.exo_names,'exact')) = cond_forecast.forecasts.controlled_exo_variables.Mean.e_a; %set controlled shocks to their values
 shock_matrix(1:5,strmatch('e_m',M_.exo_names,'exact')) = cond_forecast.forecasts.controlled_exo_variables.Mean.e_m; %set controlled shocks to their values
 
-y_simult = simult_(initial_condition_states,oo_.dr,shock_matrix,1);
+y_simult = simult_(M_,options_,initial_condition_states,oo_.dr,shock_matrix,1);
 
 if max(abs(y_simult(strmatch('k',M_.endo_names,'exact'),:)'-cond_forecast.forecasts.cond.Mean.k))>1e-8
     error('Unconditional Forecasts do not match')
diff --git a/tests/conditional_forecasts/4/fs2000_conditional_forecast_histval.mod b/tests/conditional_forecasts/4/fs2000_conditional_forecast_histval.mod
index a5c6c9dcd6188e8e4c2b59f9a2e070a192d38165..fa4912d010edb319b49e43ec36dafca5b6d4b36f 100644
--- a/tests/conditional_forecasts/4/fs2000_conditional_forecast_histval.mod
+++ b/tests/conditional_forecasts/4/fs2000_conditional_forecast_histval.mod
@@ -98,7 +98,7 @@ shock_matrix = zeros(options_cond_fcst_.periods ,M_.exo_nbr); %create shock matr
 shock_matrix(1:5,strmatch('e_a',M_.exo_names,'exact')) = cond_forecast.forecasts.controlled_exo_variables.Mean.e_a; %set controlled shocks to their values
 shock_matrix(1:5,strmatch('e_m',M_.exo_names,'exact')) = cond_forecast.forecasts.controlled_exo_variables.Mean.e_m; %set controlled shocks to their values
 
-y_simult = simult_(initial_condition_states,oo_.dr,shock_matrix,1);
+y_simult = simult_(M_,options_,initial_condition_states,oo_.dr,shock_matrix,1);
 
 if max(abs(y_simult(strmatch('k',M_.endo_names,'exact'),:)'-cond_forecast.forecasts.cond.Mean.k))>1e-8
     error('Unconditional Forecasts do not match')
diff --git a/tests/kalman_filter_smoother/compare_results_simulation/fs2000.mod b/tests/kalman_filter_smoother/compare_results_simulation/fs2000.mod
index 04a3cc7048bf46c56a812c16b70f02305acd1467..9892b9520160b18e6ffbc37beab5b2436aa21f34 100644
--- a/tests/kalman_filter_smoother/compare_results_simulation/fs2000.mod
+++ b/tests/kalman_filter_smoother/compare_results_simulation/fs2000.mod
@@ -17,7 +17,7 @@
  */
 
 /*
- * Copyright (C) 2004-2018 Dynare Team
+ * Copyright (C) 2004-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -133,7 +133,7 @@ end;
 
 dr = oo_.dr;
 iorder=1;
-y_=simult_(y0,dr,ex_,iorder);
+y_=simult_(M_,options_,y0,dr,ex_,iorder);
 
 fsdat_simul_logged;
 
diff --git a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML.mod b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML.mod
index e0a47ca36d172873b0fb90a8c4079195d4549198..fcc76dd2e2fed15c53529914026731579a30bb79 100644
--- a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML.mod
+++ b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML.mod
@@ -17,7 +17,7 @@
  */
 
 /*
- * Copyright (C) 2004-2018 Dynare Team
+ * Copyright (C) 2004-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -137,7 +137,7 @@ end;
 dr = oo_.dr;
 iorder=1;
 %run simulation 
-y_=simult_(y0,dr,ex_,iorder);
+y_=simult_(M_,options_,y0,dr,ex_,iorder);
 
 fsdat_simul_logged;
 
diff --git a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML_loglinear.mod b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML_loglinear.mod
index 473846baf8fec6889bcddbc59e63b3444369c330..13de5707551de77c052951a6ecb734cb73ecd07f 100644
--- a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML_loglinear.mod
+++ b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML_loglinear.mod
@@ -17,7 +17,7 @@
  */
 
 /*
- * Copyright (C) 2004-2018 Dynare Team
+ * Copyright (C) 2004-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -137,7 +137,7 @@ end;
 dr = oo_.dr;
 iorder=1;
 %run simulation 
-y_=simult_(y0,dr,ex_,iorder);
+y_=simult_(M_,options_,y0,dr,ex_,iorder);
 
 fsdat_simul_logged;
 
diff --git a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_loglinear.mod b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_loglinear.mod
index ff054006c3fa7370e1f0760fc704fa95713c6944..1395ddfd469acdd213b019ef6abe068e42d32f2a 100644
--- a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_loglinear.mod
+++ b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_loglinear.mod
@@ -17,7 +17,7 @@
  */
 
 /*
- * Copyright (C) 2004-2018 Dynare Team
+ * Copyright (C) 2004-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -153,7 +153,7 @@ iorder=1;
 % if options_.loglinear
 %     y0=exp(y0);
 % end
-y_=simult_(y0,dr,ex_,iorder);
+y_=simult_(M_,options_,y0,dr,ex_,iorder);
 
 fsdat_simul_logged;
 
diff --git a/tests/moments/fs2000_post_moments.mod b/tests/moments/fs2000_post_moments.mod
index dd75b5f9f7afb441b596ab1b214007d4ce4d29ac..f25a60b380d92f13349b512d322e1d8fe5c79c4a 100644
--- a/tests/moments/fs2000_post_moments.mod
+++ b/tests/moments/fs2000_post_moments.mod
@@ -130,7 +130,7 @@ par=load([M_.fname filesep 'metropolis' filesep M_.fname '_posterior_draws1']);
 
 for par_iter=1:size(par.pdraws,1)
    M_=set_parameters_locally(M_,par.pdraws{par_iter,1});
-   info=stoch_simul(var_list_);
+   [info, oo_, options_]=stoch_simul(M_, options_, oo_, var_list_);
    correlation(:,:,par_iter)=cell2mat(oo_.autocorr);
    covariance(:,:,par_iter)=oo_.var;
    conditional_variance_decomposition(:,:,:,par_iter)=oo_.conditional_variance_decomposition;
diff --git a/tests/observation_trends_and_prefiltering/Trend_no_prefilter_conditional_forecast.inc b/tests/observation_trends_and_prefiltering/Trend_no_prefilter_conditional_forecast.inc
index b2a9fca735e9fbf67dbe32a9a093c4343de303ca..c2cd41cb3630d33030489087324384e012a7ba24 100644
--- a/tests/observation_trends_and_prefiltering/Trend_no_prefilter_conditional_forecast.inc
+++ b/tests/observation_trends_and_prefiltering/Trend_no_prefilter_conditional_forecast.inc
@@ -27,7 +27,7 @@ shock_matrix = zeros(options_cond_fcst_.periods ,M_.exo_nbr); %create shock matr
 shock_matrix(1:5,strmatch('e_y',M_.exo_names,'exact')) = cond_forecast.forecasts.controlled_exo_variables.Mean.e_y; %set controlled shocks to their values
 shock_matrix(1:5,strmatch('e_p',M_.exo_names,'exact')) = cond_forecast.forecasts.controlled_exo_variables.Mean.e_p; %set controlled shocks to their values
 
-y_simult = simult_(initial_condition_states,oo_.dr,shock_matrix,1);
+y_simult = simult_(M_,options_,initial_condition_states,oo_.dr,shock_matrix,1);
 
 if max(abs(y_simult(strmatch('Y_obs',M_.endo_names,'exact'),:)'+(options_.first_obs-1+options_.nobs:options_.first_obs-1+options_.nobs+options_.forecast)'*g_y-cond_forecast.forecasts.cond.Mean.Y_obs))>1e-8
     error('Conditional Forecasts do not match')