diff --git a/matlab/simult.m b/matlab/simult.m
index 31c594cf53171da1edc106cb5957ddc61e5b980f..5c7dd9b7c9fbabffa7e0a4de2bde43b6f6cf1fb8 100644
--- a/matlab/simult.m
+++ b/matlab/simult.m
@@ -1,4 +1,4 @@
-function [y_,DynareResults] =simult(y0, dr,DynareModel,DynareOptions,DynareResults)
+function [y_out,DynareResults] =simult(y0, dr,DynareModel,DynareOptions,DynareResults)
 % Simulate a DSGE model (perturbation approach).
 
 %@info:
@@ -25,7 +25,7 @@ function [y_,DynareResults] =simult(y0, dr,DynareModel,DynareOptions,DynareResul
 %! @strong{Outputs}
 %! @sp 1
 %! @table @ @var
-%! @item y_
+%! @item y_out
 %! Matrix of doubles, simulated time series for all the endogenous variables (one per row).
 %! @item DynareResults
 %! Matlab's structure gathering the results (see @ref{oo_}).
@@ -88,6 +88,9 @@ for i=1:replic
     if replic > 1
         fwrite(fh,y_,'float64');
     end
+    if i==1
+        y_out=y_;
+    end
 end
 
 if replic > 1
diff --git a/tests/example2.mod b/tests/example2.mod
index 6022f3e8a62f42966a1568f4d0499863da34eb67..f3f7fd9407a649e60a01897c95b35b32b3cae290 100644
--- a/tests/example2.mod
+++ b/tests/example2.mod
@@ -41,3 +41,11 @@ var u = 0.009^2;
 end;
 
 stoch_simul(periods=2000, drop=200);
+
+%% test that simul_replic does not affect simulated moments
+moments_temp=oo_.var;
+set_dynare_seed('default');
+stoch_simul(periods=2000, drop=200,simul_replic=2);
+if ~isequal(moments_temp,oo_.var)
+    error('Simul_replic affects simulated moments')
+end
\ No newline at end of file