diff --git a/matlab/dyn_forecast.m b/matlab/dyn_forecast.m
index 235b4d38b39972f8fca69f18d39ff0a2a8f5642c..4343e2c8fc6e6d1511220a261064f80891ec2501 100644
--- a/matlab/dyn_forecast.m
+++ b/matlab/dyn_forecast.m
@@ -61,10 +61,10 @@ switch task
     if horizon == 0
         horizon = 5;
     end
-    if size(oo_.endo_simul,2) < maximum_lag
+    if isempty(M_.endo_histval)
         y0 = repmat(oo_.steady_state,1,maximum_lag);
     else
-        y0 = oo_.endo_simul(:,1:maximum_lag);
+        y0 = M_.endo_histval;
     end
   case 'smoother'
     horizon = options_.forecast;
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 87bd4d9acff0af61c70f83f78f50136c6aa98020..7cfb6f52878e01534fa1b9743ef977272e6c2772 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -332,6 +332,7 @@ oo_.exo_det_steady_state = [];
 oo_.exo_det_simul = [];
 
 M_.params = [];
+M_.endo_histval = [];
 
 % BVAR
 M_.bvar = [];
diff --git a/matlab/make_y_.m b/matlab/make_y_.m
index 7b6316d1ad4e8605729144b14595f355075f876e..9c1a963ed520cd63c1da7f3f62d42c686ea3c54f 100644
--- a/matlab/make_y_.m
+++ b/matlab/make_y_.m
@@ -35,13 +35,16 @@ if isempty(oo_.steady_state)
     oo_.steady_state = zeros(M_.endo_nbr,1);
 end
 
-if isempty(oo_.endo_simul)
+if isempty(M_.endo_histval)
     if isempty(ys0_)
         oo_.endo_simul = [oo_.steady_state*ones(1,M_.maximum_lag+options_.periods+M_.maximum_lead)];
     else
         oo_.endo_simul = [ys0_*ones(1,M_.maximum_lag) oo_.steady_state*ones(1,options_.periods+M_.maximum_lead)];
     end
-elseif size(oo_.endo_simul,2) < M_.maximum_lag+M_.maximum_lead+options_.periods
-        oo_.endo_simul = [oo_.endo_simul ...
-                          oo_.steady_state*ones(1,M_.maximum_lag+options_.periods+M_.maximum_lead-size(oo_.endo_simul,2),1)];
+else
+    if ~isempty(ys0_)
+        error('histval and endval cannot be used simultaneously')
+    end
+    oo_.endo_simul = [M_.endo_histval ...
+                      oo_.steady_state*ones(1,options_.periods+M_.maximum_lead)];
 end
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index ec6d0966489f20a869bcc4d215cf853d0d4dbb12..d6eb5d6257b9f85bfdb7509921b8fc81dd7fe411 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -130,13 +130,10 @@ if options_.periods > 0 && ~PI_PCL_solver
         options_ =options_old;
         return
     end
-    % Note that the first column of oo_.endo_simul is preserved by the following
-    % call to simult; this is important because stoch_simul can be followed by
-    % forecast (see ticket #157)
-    if size(oo_.endo_simul,2) == 0
+    if isempty(M_.endo_histval)
         y0 = oo_.dr.ys;
     else
-        y0 = oo_.endo_simul(:,1);
+        y0 = M_.endo_histval;
     end
     oo_.endo_simul = simult(y0,oo_.dr);
     dyn2vec;
diff --git a/matlab/stoch_simul_sparse.m b/matlab/stoch_simul_sparse.m
index 29488b9bbbc394b7c844f3779eff32d624a133fb..2136961b16f9774fb38de9e2be348b029521063e 100644
--- a/matlab/stoch_simul_sparse.m
+++ b/matlab/stoch_simul_sparse.m
@@ -82,10 +82,10 @@ elseif options_.periods ~= 0
         options_ =options_old;
         return
     end
-    if size(oo_.endo_simul,2) < maximum_lag
+    if isempty(M_.endo_histval)
         y0 = oo_.dr.ys;
     else
-        y0 = oo_.endo_simul(:,1);
+        y0 = M_.endo_histval;
     end
     oo_.endo_simul = simult(y0,oo_.dr);
     dyn2vec;
diff --git a/preprocessor/NumericalInitialization.cc b/preprocessor/NumericalInitialization.cc
index 6b08a74eb4c3292beb4dc6d74586cf344b8a9af5..7a67c56f5337c65635de2db8646b880af86fa0a9 100644
--- a/preprocessor/NumericalInitialization.cc
+++ b/preprocessor/NumericalInitialization.cc
@@ -132,8 +132,7 @@ InitValStatement::writeOutput(ostream &output, const string &basename) const
 void
 InitValStatement::writeOutputPostInit(ostream &output) const
 {
-  output << "oo_.endo_simul=[oo_.steady_state*ones(1,M_.maximum_lag)];" << endl
-         << "if M_.exo_nbr > 0;" << endl
+  output << "if M_.exo_nbr > 0;" << endl
          << "\too_.exo_simul = [ones(M_.maximum_lag,1)*oo_.exo_steady_state'];" << endl
          <<"end;" << endl
          << "if M_.exo_det_nbr > 0;" << endl
@@ -189,7 +188,7 @@ HistValStatement::writeOutput(ostream &output, const string &basename) const
   output << "%" << endl
          << "% HISTVAL instructions" << endl
          << "%" << endl
-         << "oo_.endo_simul = zeros(M_.endo_nbr,M_.maximum_lag);" << endl;
+         << "M_.endo_histval = zeros(M_.endo_nbr,M_.maximum_lag);" << endl;
 
   for (hist_values_t::const_iterator it = hist_values.begin();
        it != hist_values.end(); it++)
@@ -225,7 +224,7 @@ HistValStatement::writeOutput(ostream &output, const string &basename) const
       int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
 
       if (type == eEndogenous)
-        output << "oo_.endo_simul( " << tsid << ", M_.maximum_lag + " << lag << ") = ";
+        output << "M_.endo_histval( " << tsid << ", M_.maximum_lag + " << lag << ") = ";
       else if (type == eExogenous)
         output << "oo_.exo_simul( M_.maximum_lag + " << lag << ", " << tsid << " ) = ";
       else if (type != eExogenousDet)
diff --git a/preprocessor/NumericalInitialization.hh b/preprocessor/NumericalInitialization.hh
index b9639af910314d09849272c91b83ce6325889744..d791acace2488c71a9228312c4ae9a5c5f2f865e 100644
--- a/preprocessor/NumericalInitialization.hh
+++ b/preprocessor/NumericalInitialization.hh
@@ -71,7 +71,7 @@ public:
   InitValStatement(const init_values_t &init_values_arg,
                    const SymbolTable &symbol_table_arg);
   virtual void writeOutput(ostream &output, const string &basename) const;
-  //! Writes initializations for oo_.endo_simul, oo_.exo_simul and oo_.exo_det_simul
+  //! Writes initializations for oo_.exo_simul and oo_.exo_det_simul
   void writeOutputPostInit(ostream &output) const;
 };
 
diff --git a/tests/histval_sto.mod b/tests/histval_sto.mod
index 830e47fa1e7db02f8a19b1d0d1f6ca6883d4d56c..db4ab3c616b9e5e390dd60726b5a5e86d658baeb 100644
--- a/tests/histval_sto.mod
+++ b/tests/histval_sto.mod
@@ -49,4 +49,7 @@ a(-1) = 0.3;
 u(-1) = 0.1;
 end;
 
-stoch_simul(nograph);
+stoch_simul(nograph, periods = 200);
+
+forecast;
+