diff --git a/matlab/discretionary_policy/discretionary_policy.m b/matlab/discretionary_policy/discretionary_policy.m
index 778d1789db328591ee1aaa01ac997397a7d97cb4..b2c67403d6d2bbf87dd1f47c83b4ee540751680e 100644
--- a/matlab/discretionary_policy/discretionary_policy.m
+++ b/matlab/discretionary_policy/discretionary_policy.m
@@ -1,5 +1,5 @@
-function [info, oo_, options_] = discretionary_policy(M_, options_, oo_, var_list)
-% function [info, oo_, options_] = discretionary_policy(M_, options_, oo_, var_list)
+function [info, oo_, options_, M_] = discretionary_policy(M_, options_, oo_, var_list)
+% function [info, oo_, options_, M_] = discretionary_policy(M_, options_, oo_, var_list)
 % INPUTS
 % - M_            [structure]     Matlab's structure describing the model (M_).
 % - options_      [structure]     Matlab's structure describing the current options (options_).
@@ -10,8 +10,9 @@ function [info, oo_, options_] = discretionary_policy(M_, options_, oo_, var_lis
 % - info          [integer]       scalar or vector, error code.
 % - oo_           [structure]     Matlab's structure containing the results (oo_).
 % - options_      [structure]     Matlab's structure describing the current options (options_).
+% - M_            [structure]     Matlab's structure describing the model (M_).
 
-% Copyright (C) 2007-2019 Dynare Team
+% Copyright (C) 2007-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -33,7 +34,7 @@ M_=discretionary_policy_initialization(M_,options_);
 origorder = options_.order;
 options_.discretionary_policy = 1;
 options_.order = 1;
-[info, oo_] = stoch_simul(M_, options_, oo_, var_list);
+[info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, var_list);
 
 if ~options_.noprint
     disp_steady_state(M_,oo_)
diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m
index 6395a51b31a7382449bd22d517dad29439138840..75d429182bdb5319a527d0aa8b553aea6a028c24 100644
--- a/matlab/identification_analysis.m
+++ b/matlab/identification_analysis.m
@@ -283,7 +283,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, oo_, options_] = stoch_simul(M_, options_, oo_, options_.varobs);
+                [info, oo_, options_, M_] = 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/osr.m b/matlab/osr.m
index 25db19bdd4e52308f30535908607376d589dc131..3a0d463f2271c3f0947c9744a28205ab3149e122 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);
 
-[~, oo_, options_] = stoch_simul(M_, options_, oo_, var_list);
+[~, oo_, options_, M_] = stoch_simul(M_, options_, oo_, var_list);
diff --git a/matlab/ramsey_policy.m b/matlab/ramsey_policy.m
index 0ce386c8a3c7d0ae92d001187f3f0fadb5c932d9..46acdb85464df17891af79b057d17f892607fe24 100644
--- a/matlab/ramsey_policy.m
+++ b/matlab/ramsey_policy.m
@@ -39,7 +39,7 @@ else
     end
 end
 
-[info, oo_, options_] = stoch_simul(M_, options_, oo_, var_list);
+[info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, var_list);
 
 oo_.steady_state = oo_.dr.ys;
 
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index 4740e70323d9797d8dbb2e753ce42d2d8a34e04b..76569b3d535792545317c47a32b8faa0228ea060 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -1,4 +1,4 @@
-function [info, oo_, options_] = stoch_simul(M_, options_, oo_, var_list)
+function [info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, var_list)
 
 % Copyright (C) 2001-2020 Dynare Team
 %
diff --git a/tests/moments/fs2000_post_moments.mod b/tests/moments/fs2000_post_moments.mod
index f25a60b380d92f13349b512d322e1d8fe5c79c4a..3ea15dcef0946954f34f1a85f03927a40fa13c09 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, oo_, options_]=stoch_simul(M_, options_, oo_, var_list_);
+   [info, oo_, options_, M_]=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;