diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m
index f610a53e8764d78007216bfa8cf8833fbac61eff..1c73e474d2c09fe418199d750bc6c3e3b54ca330 100644
--- a/matlab/ep/extended_path.m
+++ b/matlab/ep/extended_path.m
@@ -1,4 +1,4 @@
-function [ts,results] = extended_path(initial_conditions,sample_size)
+function [ts,results] = extended_path(initial_conditions,sample_size, DynareOptions, DynareModel, DynareResults)
 % Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time
 % series of size T  is obtained by solving T perfect foresight models.
 %
@@ -30,36 +30,35 @@ function [ts,results] = extended_path(initial_conditions,sample_size)
 %
 % 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_
 
-ep  = options_.ep;
-options_.verbosity = ep.verbosity;
+ep  = DynareOptions.ep;
+DynareOptions.verbosity = ep.verbosity;
 verbosity = ep.verbosity+ep.debug;
 
 % Set maximum number of iterations for the deterministic solver.
-options_.simul.maxit = ep.maxit;
+DynareOptions.simul.maxit = ep.maxit;
 
 % Prepare a structure needed by the matlab implementation of the perfect foresight model solver
-pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_);
+pfm = setup_stochastic_perfect_foresight_model_solver(DynareModel,DynareOptions,DynareResults);
 
-if M_.exo_det_nbr~=0
+if DynareModel.exo_det_nbr~=0
     error('ep: Extended path does not support varexo_det.')
 end
 
-endo_nbr = M_.endo_nbr;
-exo_nbr = M_.exo_nbr;
-maximum_lag = M_.maximum_lag;
-maximum_lead = M_.maximum_lead;
+endo_nbr = DynareModel.endo_nbr;
+exo_nbr = DynareModel.exo_nbr;
+maximum_lag = DynareModel.maximum_lag;
+maximum_lead = DynareModel.maximum_lead;
 epreplic_nbr = ep.replic_nbr;
-steady_state = oo_.steady_state;
-dynatol = options_.dynatol;
+steady_state = DynareResults.steady_state;
+dynatol = DynareOptions.dynatol;
 
 % Set default initial conditions.
 if isempty(initial_conditions)
-    if isempty(M_.endo_histval)
+    if isempty(DynareModel.endo_histval)
         initial_conditions = steady_state;
     else
-        initial_conditions = M_.endo_histval;
+        initial_conditions = DynareModel.endo_histval;
     end
 end
 
@@ -68,13 +67,13 @@ end
 periods = ep.periods;
 pfm.periods = periods;
 pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
-pfm.block = options_.block;
+pfm.block = DynareOptions.block;
 
 % keep a copy of pfm.i_upd
 i_upd = pfm.i_upd;
 
 % Set the algorithm for the perfect foresight solver
-options_.stack_solve_algo = ep.stack_solve_algo;
+DynareOptions.stack_solve_algo = ep.stack_solve_algo;
 
 % Set check_stability flag
 do_not_check_stability_flag = ~ep.check_stability;
@@ -86,23 +85,23 @@ do_not_check_stability_flag = ~ep.check_stability;
 
 dr = struct();
 if ep.init
-    options_.order = 1;
-    oo_.dr=set_state_space(dr,M_,options_);
-    [dr,Info,M_,options_,oo_] = resol(0,M_,options_,oo_);
+    DynareOptions.order = 1;
+    DynareResults.dr=set_state_space(dr,DynareModel,DynareOptions);
+    [dr,Info,DynareModel,DynareOptions,DynareResults] = resol(0,DynareModel,DynareOptions,DynareResults);
 end
 
 % Do not use a minimal number of perdiods for the perfect foresight solver (with bytecode and blocks)
-options_.minimal_solving_period = 100;%options_.ep.periods;
+DynareOptions.minimal_solving_period = 100;%DynareOptions.ep.periods;
 
 % Initialize the output array.
-time_series = zeros(M_.endo_nbr,sample_size);
+time_series = zeros(DynareModel.endo_nbr,sample_size);
 
 % Set the covariance matrix of the structural innovations.
-variances = diag(M_.Sigma_e);
+variances = diag(DynareModel.Sigma_e);
 positive_var_indx = find(variances>0);
 effective_number_of_shocks = length(positive_var_indx);
 stdd = sqrt(variances(positive_var_indx));
-covariance_matrix = M_.Sigma_e(positive_var_indx,positive_var_indx);
+covariance_matrix = DynareModel.Sigma_e(positive_var_indx,positive_var_indx);
 covariance_matrix_upper_cholesky = chol(covariance_matrix);
 
 % (re)Set exo_nbr
@@ -127,12 +126,12 @@ switch ep.innovation_distribution
     shocks(:,positive_var_indx) = shocks;
   case 'calibrated'
     replic_nbr = 1;
-    shocks = zeros(sample_size,M_.exo_nbr);
-    for i = 1:length(M_.unanticipated_det_shocks)
-        k = M_.unanticipated_det_shocks(i).periods;
-        ivar = M_.unanticipated_det_shocks(i).exo_id;
-        v = M_.unanticipated_det_shocks(i).value;
-        if ~M_.unanticipated_det_shocks(i).multiplicative
+    shocks = zeros(sample_size,DynareModel.exo_nbr);
+    for i = 1:length(DynareModel.unanticipated_det_shocks)
+        k = DynareModel.unanticipated_det_shocks(i).periods;
+        ivar = DynareModel.unanticipated_det_shocks(i).exo_id;
+        v = DynareModel.unanticipated_det_shocks(i).value;
+        if ~DynareModel.unanticipated_det_shocks(i).multiplicative
             shocks(k,ivar) = v;
         else
             socks(k,ivar) = shocks(k,ivar) * v;
@@ -149,20 +148,20 @@ set(hh,'Name','EP simulations.');
 % hybrid correction
 pfm.hybrid_order = ep.stochastic.hybrid_order;
 if pfm.hybrid_order
-    oo_.dr = set_state_space(oo_.dr,M_,options_);
-    options = options_;
+    DynareResults.dr = set_state_space(DynareResults.dr,DynareModel,DynareOptions);
+    options = DynareOptions;
     options.order = pfm.hybrid_order;
-    pfm.dr = resol(0,M_,options,oo_);
+    pfm.dr = resol(0,DynareModel,options,DynareResults);
 else
     pfm.dr = [];
 end
 
 % number of nonzero derivatives
-pfm.nnzA = M_.NNZDerivatives(1);
+pfm.nnzA = DynareModel.NNZDerivatives(1);
 
 % setting up integration nodes if order > 0
 if ep.stochastic.order > 0
-    [nodes,weights,nnodes] = setup_integration_nodes(options_.ep,pfm);
+    [nodes,weights,nnodes] = setup_integration_nodes(DynareOptions.ep,pfm);
     pfm.nodes = nodes;
     pfm.weights = weights; 
     pfm.nnodes = nnodes;
@@ -175,17 +174,17 @@ end
 
 
 % set boundaries if mcp
-[lb,ub,pfm.eq_index] = get_complementarity_conditions(M_, options_.ramsey_policy);
-options_.lmmcp.lb = repmat(lb,block_nbr,1);
-options_.lmmcp.ub = repmat(ub,block_nbr,1);
+[lb,ub,pfm.eq_index] = get_complementarity_conditions(DynareModel, DynareOptions.ramsey_policy);
+DynareOptions.lmmcp.lb = repmat(lb,block_nbr,1);
+DynareOptions.lmmcp.ub = repmat(ub,block_nbr,1);
 pfm.block_nbr = block_nbr;
 
 % storage for failed draws
-oo_.ep.failures.periods = [];
-oo_.ep.failures.previous_period = cell(0);
-oo_.ep.failures.shocks = cell(0);
+DynareResults.ep.failures.periods = [];
+DynareResults.ep.failures.previous_period = cell(0);
+DynareResults.ep.failures.shocks = cell(0);
 
-oo_.exo_simul = shocks;
+DynareResults.exo_simul = shocks;
 
 % Initializes some variables.
 t  = 1;
@@ -196,7 +195,7 @@ for k = 1:replic_nbr
 end
 %make_ex_;
 exo_simul_ = zeros(maximum_lag+sample_size+maximum_lead,exo_nbr);
-exo_simul_(1:size(oo_.exo_simul,1),1:size(oo_.exo_simul,2)) = oo_.exo_simul;
+exo_simul_(1:size(DynareResults.exo_simul,1),1:size(DynareResults.exo_simul,2)) = DynareResults.exo_simul;
 % Main loop.
 while (t <= sample_size)
     if ~mod(t,10)
@@ -207,21 +206,21 @@ while (t <= sample_size)
    
     if replic_nbr > 1 && ep.parallel_1
         parfor k = 1:replic_nbr
-            exo_simul = repmat(oo_.exo_steady_state',periods+2,1);
+            exo_simul = repmat(DynareResults.exo_steady_state',periods+2,1);
             %            exo_simul(1:sample_size+3-t,:) = exo_simul_(t:end,:);
-            exo_simul(2,:) = exo_simul_(M_.maximum_lag+t,:) + ...
+            exo_simul(2,:) = exo_simul_(DynareModel.maximum_lag+t,:) + ...
                 shocks((t-2)*replic_nbr+k,:);
             initial_conditions = results{k}(:,t-1);
             [results{k}(:,t), info_convergence] = extended_path_core(ep.periods,endo_nbr,exo_nbr,positive_var_indx, ...
                                                               exo_simul,ep.init,initial_conditions,...
                                                               maximum_lag,maximum_lead,steady_state, ...
                                                               ep.verbosity,bytecode_flag,ep.stochastic.order,...
-                                                              M_.params,pfm,ep.stochastic.algo,ep.solve_algo,ep.stack_solve_algo,...
-                                                              options_.lmmcp,options_,oo_);
+                                                              DynareModel.params,pfm,ep.stochastic.algo,ep.solve_algo,ep.stack_solve_algo,...
+                                                              DynareOptions.lmmcp,DynareOptions,DynareResults);
         end
     else
         for k = 1:replic_nbr
-            exo_simul = repmat(oo_.exo_steady_state',periods+maximum_lag+ ...
+            exo_simul = repmat(DynareResults.exo_steady_state',periods+maximum_lag+ ...
                             maximum_lead,1);
             %            exo_simul(1:sample_size+maximum_lag+maximum_lead-t+1,:) = ...
             %                exo_simul_(t:end,:);
@@ -232,8 +231,8 @@ while (t <= sample_size)
                                                               exo_simul,ep.init,initial_conditions,...
                                                               maximum_lag,maximum_lead,steady_state, ...
                                                               ep.verbosity,bytecode_flag,ep.stochastic.order,...
-                                                              M_,pfm,ep.stochastic.algo,ep.solve_algo,ep.stack_solve_algo,...
-                                                              options_.lmmcp,options_,oo_);
+                                                              DynareModel,pfm,ep.stochastic.algo,ep.solve_algo,ep.stack_solve_algo,...
+                                                              DynareOptions.lmmcp,DynareOptions,DynareResults);
         end
     end
     if verbosity
@@ -247,25 +246,25 @@ end% (while) loop over t
 
 dyn_waitbar_close(hh);
 
-if isnan(options_.initial_period)
+if isnan(DynareOptions.initial_period)
     initial_period = dates(1,1);
 else
-    initial_period = options_.initial_period;
+    initial_period = DynareOptions.initial_period;
 end
 if nargout
     if ~isnan(results{1})
         ts = dseries(transpose([results{1}]), ...
-                     initial_period,cellstr(M_.endo_names));
+                     initial_period,cellstr(DynareModel.endo_names));
     else
         ts = NaN;
     end
 else
     if ~isnan(results{1})
-        oo_.endo_simul = results{1};
+        DynareResults.endo_simul = results{1};
         ts = dseries(transpose(results{1}),initial_period, ...
-                     cellstr(M_.endo_names));
+                     cellstr(DynareModel.endo_names));
     else
-        oo_.endo_simul = NaN;
+        DynareResults.endo_simul = NaN;
         ts = NaN;
     end
 end
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index 16cfde025e24adfc8a2905d1e8cefee91c0aee6f..24c3f12d3fa7af0bb7c444971f14d8acc196ea41 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -3200,7 +3200,7 @@ ExtendedPathStatement::writeOutput(ostream &output, const string &basename, bool
       output << "options_." << it->first << " = " << it->second << ";" << endl;
 
   output << "extended_path([], " << options_list.num_options.find("periods")->second
-         << ");" << endl;
+         << ", options_, M_, oo_);" << endl;
 }
 
 ModelDiagnosticsStatement::ModelDiagnosticsStatement()
diff --git a/tests/ep/ar.mod b/tests/ep/ar.mod
index 4945a883ebc861444ce5b59f32093e24d4a351b9..4cc34cf4ad3a2e99808cd9747e7208866e492334 100644
--- a/tests/ep/ar.mod
+++ b/tests/ep/ar.mod
@@ -36,7 +36,7 @@ options_.ep.stochastic.order = 0;
 options_.ep.stochastic.nodes = 0;
 options_.console_mode = 0;
 
-ts = extended_path([],10);
+ts = extended_path([], 10, options_, M_, oo_);
 
 options_.ep.verbosity = 0;
 options_.ep.stochastic.order = 1;
@@ -44,7 +44,7 @@ options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature';
 options_.ep.stochastic.nodes = 3;
 options_.console_mode = 0;
 
-sts = extended_path([],10);
+sts = extended_path([], 10, options_, M_, oo_);
 
 if max(max(abs(ts-sts)))>pi*options_.dynatol.x
    disp('Stochastic Extended Path:: Something is wrong here (potential bug in extended_path.m)!!!')
diff --git a/tests/ep/burnside.mod b/tests/ep/burnside.mod
index 7253d2f0b8cd7658d214336da3499315d39300a7..cd863b8215134f5bfa6ebd23a579eb42e6895058 100644
--- a/tests/ep/burnside.mod
+++ b/tests/ep/burnside.mod
@@ -50,15 +50,15 @@ options_.ep.stochastic.nodes = 2;
 options_.console_mode = 0;
 
 set_dynare_seed('default');
-ts = extended_path([],5000);
+ts = extended_path([], 5000, options_, M_, oo_);
 
 options_.ep.stochastic.order = 2;
 options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature';
 set_dynare_seed('default');
-ts1_4 = extended_path([],5000);
+ts1_4 = extended_path([], 5000, options_, M_, oo_);
 
 set_dynare_seed('default');
-ytrue=exact_solution(M_,oo_,800);
+ytrue=exact_solution(M_,oo_, 800);
 
 disp('True mean and standard deviation')
 disp(mean(ytrue(101:end)))
diff --git a/tests/ep/linearmodel.mod b/tests/ep/linearmodel.mod
index 1de421589a47637802bf1df584832fbea17aa9f6..a36902dd6c4ddc6195febe6daec2b729526b4280 100644
--- a/tests/ep/linearmodel.mod
+++ b/tests/ep/linearmodel.mod
@@ -33,13 +33,13 @@ options_.ep.order = 0;
 options_.ep.nnodes = 0;
 options_.console_mode = 0;
 
-ts = extended_path([],10);
+ts = extended_path([], 10, options_, M_, oo_);
 
 options_.ep.stochastic.status = 1;
 options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature';
 options_.ep.order = 1;
 options_.ep.nnodes = 3;
-sts = extended_path([],10);
+sts = extended_path([], 10, options_, M_, oo_);
 
 if max(max(abs(ts.data-sts.data))) > 1e-12
    error('extended path algorithm fails in ./tests/ep/linearmodel.mod')
diff --git a/tests/ep/rbc.mod b/tests/ep/rbc.mod
index 5c608b6a0ac1fdc264bdc1a0230289a914df2309..4c2e3dc475b5b85984a9cc27075e20370bc47980 100644
--- a/tests/ep/rbc.mod
+++ b/tests/ep/rbc.mod
@@ -75,19 +75,19 @@ steady(nocheck);
 options_.ep.verbosity = 0;
 
 options_.ep.stochastic.order = 0;
-ts0 = extended_path([],10);
+ts0 = extended_path([], 10, options_, M_, oo_);
 
 options_.ep.stochastic.order = 1;
 options_.ep.stochastic.nodes = 3;
 options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature';
-ts1_3 = extended_path([],10);
+ts1_3 = extended_path([], 10, options_, M_, oo_);
 
 options_.ep.stochastic.nodes = 5;
-ts1_5 = extended_path([],10);
+ts1_5 = extended_path([], 10, options_, M_, oo_);
 
 options_.ep.stochastic.order = 2;
 options_.ep.stochastic.nodes = 3;
-ts2_3 = extended_path([],10);
+ts2_3 = extended_path([], 10, options_, M_, oo_);
 
 options_.ep.stochastic.nodes = 5;
-ts2_5 = extended_path([],10);
+ts2_5 = extended_path([], 10, options_, M_, oo_);
diff --git a/tests/ep/rbcii.mod b/tests/ep/rbcii.mod
index 18f7a9bb72f55c5c4cc50d616410b827be47d921..2fa42bffa1daa98658afc16b7cf25e83fe897d6f 100644
--- a/tests/ep/rbcii.mod
+++ b/tests/ep/rbcii.mod
@@ -72,12 +72,12 @@ copyfile('rbcii_steady_state.m','rbcii_steadystate2.m');
     options_.ep.stochastic.nodes = 2;
     options_.console_mode = 0;
 
-    ts = extended_path([],20);
+    ts = extended_path([], 20, options_, M_, oo_);
 
     options_.ep.stochastic.order = 1;
     options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature';
 //    profile on
-    ts1_4 = extended_path([],20);
+    ts1_4 = extended_path([], 20, options_, M_, oo_);
 //    profile off
 //    profile viewer
 @#else