diff --git a/mex/sources/ms-sbvar/mex_ms_forecast.cc b/mex/sources/ms-sbvar/mex_ms_forecast.cc
index 678be6b95..08d5eabbe 100644
--- a/mex/sources/ms-sbvar/mex_ms_forecast.cc
+++ b/mex/sources/ms-sbvar/mex_ms_forecast.cc
@@ -115,6 +115,7 @@ mexFunction(int nlhs, mxArray *plhs[],
         for (s = 0; s < nstates; s++)
           {
             printf("Constructing percentiles for forecast - regime %d\n", s);
+            if (options->simulation_file) rewind(options->simulation_file);
             if (histogram = forecast_percentile_regime(options->shocks, options->simulation_file, options->thin, s, T, options->horizon, model, type))
               {
                 for (k = 0; k < options->num_percentiles; k++)
diff --git a/mex/sources/ms-sbvar/mex_ms_irf.cc b/mex/sources/ms-sbvar/mex_ms_irf.cc
index ae607dd57..fe0131562 100644
--- a/mex/sources/ms-sbvar/mex_ms_irf.cc
+++ b/mex/sources/ms-sbvar/mex_ms_irf.cc
@@ -114,6 +114,7 @@ mexFunction(int nlhs, mxArray *plhs[],
         for (s = 0; s < nstates; s++)
           {
             printf("Constructing percentiles for impulse responses - regime %d\n", s);
+            if (options->simulation_file) rewind(options->simulation_file);
             if (histogram = impulse_response_percentile_regime(options->simulation_file, options->thin, s, options->horizon, model, type))
               {
                 for (k = 0; k < options->num_percentiles; k++)
diff --git a/mex/sources/ms-sbvar/mex_ms_variance_decomposition.cc b/mex/sources/ms-sbvar/mex_ms_variance_decomposition.cc
index fdf0df038..9d8e813a6 100644
--- a/mex/sources/ms-sbvar/mex_ms_variance_decomposition.cc
+++ b/mex/sources/ms-sbvar/mex_ms_variance_decomposition.cc
@@ -130,6 +130,7 @@ mexFunction(int nlhs, mxArray *plhs[],
                 }
               else
                 {
+                  if (options->simulation_file) rewind(options->simulation_file);
                   if (histogram = variance_decomposition_percentiles_regime(options->simulation_file, options->thin, s, options->horizon, model, type))
                     {
                       tvd = CreateMatrix(options->horizon, nvars*nvars);