diff --git a/doc/dynare.texi b/doc/dynare.texi
index e3049ec56b704df72e977600fd0eff2759a0267d..944c952568e254cea4f46287f77abb9bd165efa9 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -5933,10 +5933,6 @@ while data files are contained in @code{<output_file_tag/Forecast>}.
 @item data_obs_nbr = @var{INTEGER}
 The number of data points included in the output. Default: @code{0}
 
-@item no_error_bands
-@anchor{no_error_bands} Do not output error bands. Default: @code{off}
-(@i{i.e.} output error bands)
-
 @item error_band_percentiles = [@var{DOUBLE1} @dots{}]
 @xref{error_band_percentiles}.
 
@@ -5952,6 +5948,16 @@ The number of data points included in the output. Default: @code{0}
 @item free_parameters = @var{NUMERICAL_VECTOR}
 @xref{free_parameters}.
 
+@item parameter_uncertainty
+@xref{parameter_uncertainty}.
+
+@item regime = @var{INTEGER}
+@xref{regime}.
+
+@item regimes
+
+@xref{regimes}.
+
 @item median
 
 @xref{median}.
@@ -5988,7 +5994,8 @@ are contained in @code{<output_file_tag/Variance_Decomposition>}.
 @xref{filtered_probabilities}.
 
 @item no_error_bands
-@xref{no_error_bands}.
+Do not output error bands. Default: @code{off}
+(@i{i.e.} output error bands)
 
 @item error_band_percentiles = [@var{DOUBLE1} @dots{}]
 @xref{error_band_percentiles}.
diff --git a/matlab/ms-sbvar/initialize_ms_sbvar_options.m b/matlab/ms-sbvar/initialize_ms_sbvar_options.m
index 3937641e8dce8dae9869758bc8844a242e8234d0..69da2963ad27f02336867f2dd46a06e1914c0d69 100644
--- a/matlab/ms-sbvar/initialize_ms_sbvar_options.m
+++ b/matlab/ms-sbvar/initialize_ms_sbvar_options.m
@@ -108,6 +108,7 @@ options_.ms.median = 0;
 options_.ms.regime = 0;
 options_.ms.regimes = 0;
 % forecast
-options_.ms.error_bands = 1;
 options_.ms.forecast_data_obs = 0;
+% variance decomposition
+options_.ms.error_bands = 1;
 end
diff --git a/matlab/ms-sbvar/ms_forecast.m b/matlab/ms-sbvar/ms_forecast.m
index a8788ba7482d23e718c4d673a0c8f5ff6b54d511..e267319b9379ac1b252f226e72471285346c6543 100644
--- a/matlab/ms-sbvar/ms_forecast.m
+++ b/matlab/ms-sbvar/ms_forecast.m
@@ -14,7 +14,7 @@ function [options_, oo_]=ms_forecast(M_, options_, oo_)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2011 Dynare Team
+% Copyright (C) 2011-2012 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -34,46 +34,88 @@ function [options_, oo_]=ms_forecast(M_, options_, oo_)
 disp('MS-SBVAR Forecasts');
 options_ = set_file_tags(options_);
 [options_, oo_] = set_ms_estimation_file(options_.ms.file_tag, options_, oo_);
-options_ = set_ms_simulation_file(options_);
-clean_files_for_second_type_of_mex(M_, options_, 'forecast')
+clean_ms_forecast_files(options_.ms.output_file_tag);
 forecastdir = [options_.ms.output_file_tag filesep 'Forecast'];
 create_dir(forecastdir);
 
-opt = { ...
-    {'file_tag', options_.ms.file_tag}, ...
-    {'seed', options_.DynareRandomStreams.seed}, ...
-    {'horizon', options_.ms.horizon}, ...
-    {'number_observations', options_.ms.forecast_data_obs}, ...
-    {'error_bands', options_.ms.error_bands}, ...
-    {'percentiles', options_.ms.percentiles}, ...
-    {'thin', options_.ms.thinning_factor}
-    };
+% setup command line options
+opt = ['-forecast -nodate -seed ' num2str(options_.DynareRandomStreams.seed)];
+opt = [opt ' -do ' forecastdir];
+opt = [opt ' -ft ' options_.ms.file_tag];
+opt = [opt ' -fto ' options_.ms.output_file_tag];
+opt = [opt ' -horizon ' num2str(options_.ms.horizon)];
+opt = [opt ' -thin ' num2str(options_.ms.thinning_factor)];
+opt = [opt ' -data ' num2str(options_.ms.forecast_data_obs)];
 
+if options_.ms.regimes
+    opt = [opt ' -regimes'];
+elseif options_.ms.regime
+    % regime-1 since regime is 0-indexed in C but 1-indexed in Matlab
+    opt = [opt ' -regime ' num2str(options_.ms.regime-1)];
+end
+
+if options_.ms.parameter_uncertainty
+    options_ = set_ms_simulation_file(options_);
+    opt = [opt ' -parameter_uncertainty'];
+    opt = [opt ' -shocks_per_parameter ' num2str(options_.ms.shocks_per_parameter)];
+else
+    opt = [opt ' -shocks_per_parameter ' num2str(options_.ms.shock_draws)];
+end
+
+percentiles_size = 0;
 if options_.ms.median
-    opt = [opt(:)' {{'median'}}];
+    percentiles_size = 1;
+    opt = [opt ' -percentiles ' num2str(percentiles_size) ' 0.5'];
+else
+    percentiles_size = size(options_.ms.percentiles,2);
+    opt = [opt ' -percentiles ' num2str(percentiles_size)];
+    for i=1:size(options_.ms.percentiles,2)
+        opt = [opt ' ' num2str(options_.ms.percentiles(i))];
+    end
 end
 
-[err, forecast] = mex_ms_forecast([opt(:)', {{'free_parameters',oo_.ms.maxparams}, ...
-    {'shocks_per_parameter', options_.ms.shock_draws}}]);
-mexErrCheck('mex_ms_forecast ergodic ', err);
-plot_ms_forecast(M_,options_,forecast,'Forecast',options_.graph_save_formats,options_.TeX);
+% forecast
+[err] = ms_sbvar_command_line(opt);
+mexErrCheck('ms_forecast',err);
 
-[err, regime_forecast] = mex_ms_forecast([opt(:)', {{'free_parameters',oo_.ms.maxparams}, ...
-    {'shocks_per_parameter', options_.ms.shock_draws}, {'regimes'}}]);
-mexErrCheck('mex_ms_forecast ergodic regimes', err);
-save([forecastdir filesep 'ergodic_forecast.mat'], 'forecast', 'regime_forecast');
+% Plot Forecasts
+if options_.ms.regimes
+    n_chains = length(options_.ms.ms_chain);
+    n_states=1;
+    for i_chain=1:n_chains
+        n_states = n_states*length(options_.ms.ms_chain(i_chain).regime);
+    end
 
-if exist(options_.ms.mh_file,'file') > 0
-    [err, forecast] = mex_ms_forecast([opt(:)', {{'free_parameters',oo_.ms.maxparams}, ...
-        {'shocks_per_parameter', options_.ms.shocks_per_parameter}, ...
-        {'simulation_file', options_.ms.mh_file}, {'parameter_uncertainty'}}]);
-    mexErrCheck('mex_ms_forecast bayesian ', err);
-    plot_ms_forecast(M_,options_,forecast,'Forecast w/ Parameter Uncertainty',options_.graph_save_formats,options_.TeX);
+    for state_i=1:n_states
+        forecast_data = load([forecastdir filesep 'forecasts_percentiles_regime_' ...
+            num2str(state_i-1) '_' options_.ms.output_file_tag ...
+            '.out'], '-ascii');
+        forecast_data = reshape_ascii_forecast_data(M_.endo_nbr, ...
+            percentiles_size, options_.ms.horizon, forecast_data);
+        save([forecastdir filesep 'forecast_state_' num2str(state_i-1)], ...
+            'forecast_data');
+        plot_ms_forecast(M_, options_, forecast_data, ...
+            ['Forecast, Regimes' num2str(state_i)], ...
+            options_.graph_save_formats, options_.TeX);
+    end
+else
+    if options_.ms.regime
+        forecast_data = load([forecastdir filesep 'forecasts_percentiles_regime_' ...
+            num2str(options_.ms.regime-1) '_' options_.ms.output_file_tag ...
+            '.out'], '-ascii');
+        forecast_title = ['Forecast, State ' num2str(options_.ms.regime)];
+        save_filename = ['forecast_regime_' num2str(options_.ms.regime-1)];
+    else
+        forecast_data = load([forecastdir filesep 'forecasts_percentiles_' ...
+            options_.ms.output_file_tag '.out'], '-ascii');
+        forecast_title = 'Forecast';
+        save_filename = 'forecast';
+    end
 
-    [err, regime_forecast] = mex_ms_forecast([opt(:)', {{'free_parameters',oo_.ms.maxparams}, ...
-        {'shocks_per_parameter', options_.ms.shocks_per_parameter}, ...
-        {'simulation_file', options_.ms.mh_file}, {'parameter_uncertainty','regimes'}}]);
-    mexErrCheck('mex_ms_forecast bayesian regimes ', err);
-    save([forecastdir filesep 'bayesian_forecast.mat'], 'forecast', 'regime_forecast');
+    forecast_data = reshape_ascii_forecast_data(M_.endo_nbr, ...
+        percentiles_size, options_.ms.horizon, forecast_data);
+    save([forecastdir filesep save_filename], 'forecast_data');
+    plot_ms_forecast(M_, options_, forecast_data, forecast_title, ...
+        options_.graph_save_formats, options_.TeX);
 end
 end
diff --git a/matlab/ms-sbvar/reshape_ascii_forecast_data.m b/matlab/ms-sbvar/reshape_ascii_forecast_data.m
new file mode 100644
index 0000000000000000000000000000000000000000..6c2c06ff40b7cd6977ff4705eebfba0df2657347
--- /dev/null
+++ b/matlab/ms-sbvar/reshape_ascii_forecast_data.m
@@ -0,0 +1,44 @@
+function forecast_data=reshape_ascii_forecast_data(endo_nbr, psize, horizon, ascii_data)
+% function forecast_data=reshape_ascii_forecast_data(endo_nbr, psize, horizon, ascii_data)
+%
+% INPUTS
+%    endo_nbr:    number of endogenous
+%    psize:       number of percentiles
+%    horizon:     forecast horizon
+%    ascii_data:  data from .out file created by Dan's C code
+%
+% OUTPUTS
+%    forecast_data:    new 3-d array holding data with error bands
+%
+% SPECIAL REQUIREMENTS
+%    none
+
+% Copyright (C) 2011-2012 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+if psize <= 1
+    forecast_data = ascii_data;
+    return
+end
+
+forecast_data = zeros(psize, horizon, endo_nbr);
+for i=1:endo_nbr
+    for j=1:psize
+        forecast_data(j,:,i) = ascii_data(1+horizon*(j-1):horizon*j,i)';
+    end
+end
+end
diff --git a/mex/build/ms_sbvar.am b/mex/build/ms_sbvar.am
index 1f703bbed7ff4e7ac022bd3163c20d3e85fb3e36..6270789b224366148fd7a810f891fcee55c31b77 100644
--- a/mex/build/ms_sbvar.am
+++ b/mex/build/ms_sbvar.am
@@ -1,4 +1,4 @@
-noinst_PROGRAMS = ms_sbvar_create_init_file ms_sbvar_command_line mex_ms_convert_free_parameters mex_ms_forecast mex_ms_variance_decomposition
+noinst_PROGRAMS = ms_sbvar_create_init_file ms_sbvar_command_line mex_ms_convert_free_parameters mex_ms_variance_decomposition
 
 DWSWITCHDIR = $(top_srcdir)/../../../ms-sbvar/switch_dw
 DWUTILITIESDIR = $(top_srcdir)/../../../ms-sbvar/utilities_dw
@@ -62,9 +62,9 @@ nodist_ms_sbvar_command_line_SOURCES = \
 	$(DWSWITCHDIR)/state_space/sbvar/dw_csminwel.c \
 	$(DWUTILS) \
 	$(MSMEXSRC)
-	
+
 nodist_mex_ms_convert_free_parameters_SOURCES = \
-  $(DWSWITCHDIR)/switching/dw_switch.c \
+	$(DWSWITCHDIR)/switching/dw_switch.c \
 	$(DWSWITCHDIR)/switching/dw_switchio.c \
 	$(DWSWITCHDIR)/switching/dw_dirichlet_restrictions.c \
 	$(DWSWITCHDIR)/switching/dw_metropolis_theta.c \
@@ -73,51 +73,28 @@ nodist_mex_ms_convert_free_parameters_SOURCES = \
 	$(DWSWITCHDIR)/state_space/sbvar/VARio_matlab.c \
 	$(MSMEXSRCDIR)/mex_ms_convert_free_parameters.cc \
 	$(MSMEXSRCDIR)/mex_ms_sbvar.cc \
-  $(MSMEXSRCDIR)/modify_for_mex.cc \
+	$(MSMEXSRCDIR)/modify_for_mex.cc \
 	$(DWUTILS)
 
-nodist_mex_ms_forecast_SOURCES = \
-  $(DWSWITCHDIR)/switching/dw_switch.c \
-  $(DWSWITCHDIR)/switching/dw_switchio.c \
-  $(DWSWITCHDIR)/switching/dw_dirichlet_restrictions.c \
-  $(DWSWITCHDIR)/switching/dw_metropolis_theta.c \
-  $(DWSWITCHDIR)/switching/dw_switch_opt.c \
-  $(DWSWITCHDIR)/switching/dw_mdd_switch.c \
-  $(DWSWITCHDIR)/state_space/sbvar/VARbase.c \
-  $(DWSWITCHDIR)/state_space/sbvar/VARio.c \
-  $(DWSWITCHDIR)/state_space/sbvar/dw_sbvar_command_line.c \
-  $(DWSWITCHDIR)/state_space/sbvar/sbvar_estimate.c \
-  $(DWSWITCHDIR)/state_space/sbvar/sbvar_simulate.c \
-  $(DWSWITCHDIR)/state_space/sbvar/sbvar_probabilities.c \
-  $(DWSWITCHDIR)/state_space/sbvar/sbvar_mdd.c \
-  $(DWSWITCHDIR)/state_space/sbvar/sbvar_forecast.c \
-  $(DWSWITCHDIR)/state_space/sbvar/sbvar_variance_decomposition.c \
-  $(DWSWITCHDIR)/state_space/sbvar/sbvar_impulse_responses.c \
-  $(DWSWITCHDIR)/state_space/sbvar/dw_csminwel.c \
-  $(MSMEXSRCDIR)/mex_ms_sbvar.cc \
-  $(MSMEXSRCDIR)/modify_for_mex.cc \
-	$(MSMEXSRCDIR)/mex_ms_forecast.cc \
-	$(DWUTILS)
-	
 nodist_mex_ms_variance_decomposition_SOURCES = \
-  $(DWSWITCHDIR)/switching/dw_switch.c \
-  $(DWSWITCHDIR)/switching/dw_switchio.c \
-  $(DWSWITCHDIR)/switching/dw_dirichlet_restrictions.c \
-  $(DWSWITCHDIR)/switching/dw_metropolis_theta.c \
-  $(DWSWITCHDIR)/switching/dw_switch_opt.c \
-  $(DWSWITCHDIR)/switching/dw_mdd_switch.c \
-  $(DWSWITCHDIR)/state_space/sbvar/VARbase.c \
-  $(DWSWITCHDIR)/state_space/sbvar/VARio.c \
-  $(DWSWITCHDIR)/state_space/sbvar/dw_sbvar_command_line.c \
-  $(DWSWITCHDIR)/state_space/sbvar/sbvar_estimate.c \
-  $(DWSWITCHDIR)/state_space/sbvar/sbvar_simulate.c \
-  $(DWSWITCHDIR)/state_space/sbvar/sbvar_probabilities.c \
-  $(DWSWITCHDIR)/state_space/sbvar/sbvar_mdd.c \
-  $(DWSWITCHDIR)/state_space/sbvar/sbvar_forecast.c \
-  $(DWSWITCHDIR)/state_space/sbvar/sbvar_variance_decomposition.c \
-  $(DWSWITCHDIR)/state_space/sbvar/sbvar_impulse_responses.c \
-  $(DWSWITCHDIR)/state_space/sbvar/dw_csminwel.c \
-  $(MSMEXSRCDIR)/mex_ms_sbvar.cc \
-  $(MSMEXSRCDIR)/modify_for_mex.cc \
+	$(DWSWITCHDIR)/switching/dw_switch.c \
+	$(DWSWITCHDIR)/switching/dw_switchio.c \
+	$(DWSWITCHDIR)/switching/dw_dirichlet_restrictions.c \
+	$(DWSWITCHDIR)/switching/dw_metropolis_theta.c \
+	$(DWSWITCHDIR)/switching/dw_switch_opt.c \
+	$(DWSWITCHDIR)/switching/dw_mdd_switch.c \
+	$(DWSWITCHDIR)/state_space/sbvar/VARbase.c \
+	$(DWSWITCHDIR)/state_space/sbvar/VARio.c \
+	$(DWSWITCHDIR)/state_space/sbvar/dw_sbvar_command_line.c \
+	$(DWSWITCHDIR)/state_space/sbvar/sbvar_estimate.c \
+	$(DWSWITCHDIR)/state_space/sbvar/sbvar_simulate.c \
+	$(DWSWITCHDIR)/state_space/sbvar/sbvar_probabilities.c \
+	$(DWSWITCHDIR)/state_space/sbvar/sbvar_mdd.c \
+	$(DWSWITCHDIR)/state_space/sbvar/sbvar_forecast.c \
+	$(DWSWITCHDIR)/state_space/sbvar/sbvar_variance_decomposition.c \
+	$(DWSWITCHDIR)/state_space/sbvar/sbvar_impulse_responses.c \
+	$(DWSWITCHDIR)/state_space/sbvar/dw_csminwel.c \
+	$(MSMEXSRCDIR)/mex_ms_sbvar.cc \
+	$(MSMEXSRCDIR)/modify_for_mex.cc \
 	$(MSMEXSRCDIR)/mex_ms_variance_decomposition.cc \
 	$(DWUTILS)
diff --git a/mex/sources/ms-sbvar/mex_ms_forecast.cc b/mex/sources/ms-sbvar/mex_ms_forecast.cc
deleted file mode 100644
index 5a27ca3a4394a06cc25fd2e690bf40f6abba8f00..0000000000000000000000000000000000000000
--- a/mex/sources/ms-sbvar/mex_ms_forecast.cc
+++ /dev/null
@@ -1,187 +0,0 @@
-/*
- * Copyright (C) 2011 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
-#include "dynmex.h"
-#include "mex_ms_sbvar.h"
-
-#include <iostream>
-
-extern "C"  {
-#include "modify_for_mex.h"
-#include "switch.h"
-#include "switchio.h"
-#include "VARio.h"
-#include "dw_histogram.h"
-#include "sbvar_forecast.h"
-}
-
-void
-mexFunction(int nlhs, mxArray *plhs[],
-            int nrhs, const mxArray *prhs[])
-{
-  double *out_buf;
-  int i, j, k, s, nfree, nstates, nvars, npre, T;
-
-  TStateModel *model;
-  SbvarOption *options = (SbvarOption *) NULL;
-  T_VAR_Parameters *p = (T_VAR_Parameters *) NULL;
-  TMatrixHistogram *histogram = (TMatrixHistogram *) NULL;
-  TMatrix forecast;
-
-  int type = F_FREE;
-
-  // Check the left hand and right hand side arguments to make sure they conform
-  if (nrhs != 1)
-    DYN_MEX_FUNC_ERR_MSG_TXT("ms_forecast takes one cell array as an input argument.");
-  if (nlhs != 2)
-    DYN_MEX_FUNC_ERR_MSG_TXT("You must specify two output arguments.");
-
-  model = initialize_model_and_options(&options, prhs, &nstates, &nvars, &npre, &nfree);
-  if (model == NULL || options == NULL)
-    DYN_MEX_FUNC_ERR_MSG_TXT("There was a problem initializing the model, can not continue");
-
-  p = (T_VAR_Parameters *) (model->theta);
-  T = p->nobs;
-
-  // Check to make sure that there is a simulation file present if we are
-  // using parameter uncertainty
-  if (options->parameter_uncertainty && options->simulation_file == NULL)
-    DYN_MEX_FUNC_ERR_MSG_TXT("Paramter Uncertainty Was Specified but the simulation file was not found, please specify the simulation file to use with: 'simulation_file',<filename>");
-  if (!options->parameter_uncertainty)
-    options->simulation_file = (FILE *) NULL;
-
-  // Allocate the output matrix
-  if (options->regimes)
-    if (options->num_percentiles > 1)
-      {
-        // regimes x percentile x horizon x nvar
-        mwSize dims[4] = {nstates, options->num_percentiles, options->horizon, nvars};
-        plhs[1] = mxCreateNumericArray(4, dims, mxDOUBLE_CLASS, mxREAL);
-        out_buf = mxGetPr(plhs[1]);
-      }
-    else
-      {
-        // regimes x horizon x (nvar*nvar)
-        mwSize dims[3] = {nstates, options->horizon, nvars};
-        plhs[1] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
-        out_buf = mxGetPr(plhs[1]);
-      }
-  else
-  if (options->num_percentiles > 1)
-    {
-      // percentile x horizon x (nvar*nvar)
-      mwSize dims[3] = {options->num_percentiles, options->horizon, nvars};
-      plhs[1] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
-      out_buf = mxGetPr(plhs[1]);
-    }
-  else
-    {
-      // horizon x (nvar*nvar)
-      mwSize dims[2] = {options->horizon, nvars};
-      plhs[1] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL);
-      out_buf = mxGetPr(plhs[1]);
-    }
-
-  // hard coding dates for now to be off
-  int dates = 0;
-  try
-    {
-      if (options->regimes)
-        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++)
-                  if (forecast = ForecastHistogramToPercentile((TMatrix) NULL, options->number_observations, options->horizon, T, options->percentiles[k], histogram, model, dates))
-                    {
-                      // the case where we are only dealing with one percentile so we can reduce dimensionality of output argument
-                      if (options->num_percentiles == 1)
-                        for (i = 0; i < options->horizon; i++)
-                          for (j = 0; j < nvars; j++)
-                            out_buf[(s) + ((i) + (j)*options->horizon)*nstates] = ElementM(forecast, i, j);
-                      else
-                        for (i = 0; i < options->horizon; i++)
-                          for (j = 0; j < nvars; j++)
-                            out_buf[(s) + ((k) + ((i) + (j)*options->horizon)*options->num_percentiles)*nstates] = ElementM(forecast, i, j);
-                      mxFree(forecast);
-                    }
-                mxFree(histogram);
-              }
-          }
-
-      else if (options->regime >= 0)
-        {
-          s = options->regime;
-          printf("Constructing percentiles for forecast - regime %d\n", s);
-          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++)
-                if (forecast = ForecastHistogramToPercentile((TMatrix) NULL, options->number_observations, options->horizon, T, options->percentiles[k], histogram, model, dates))
-                  {
-                    // the case where we are only dealing with one percentile so we can reduce dimensionality of output argument
-                    if (options->num_percentiles == 1)
-                      for (i = 0; i < options->horizon; i++)
-                        for (j = 0; j < nvars; j++)
-                          out_buf[(i) + (j)*options->horizon] = ElementM(forecast, i, j);
-                    else
-                      for (i = 0; i < options->horizon; i++)
-                        for (j = 0; j < nvars; j++)
-                          out_buf[(k) + ((i) + (j)*options->horizon)*options->num_percentiles] = ElementM(forecast, i, j);
-                    mxFree(forecast);
-                  }
-              mxFree(histogram);
-            }
-        }
-      else
-        {
-          printf("Constructing percentiles for forecasts - %d draws of shocks/regimes per posterior value\n", options->shocks);
-          if (histogram = forecast_percentile(options->shocks, options->simulation_file, options->thin, T, options->horizon, model, type))
-            {
-              for (k = 0; k < options->num_percentiles; k++)
-                if (forecast = ForecastHistogramToPercentile((TMatrix) NULL, options->number_observations, options->horizon, T, options->percentiles[k], histogram, model, dates))
-                  {
-                    // the case where we are only dealing with one percentile so we can reduce dimensionality of output argument
-                    if (options->num_percentiles == 1)
-                      for (i = 0; i < options->horizon; i++)
-                        for (j = 0; j < nvars; j++)
-                          out_buf[(i) + (j)*options->horizon] = ElementM(forecast, i, j);
-                    else
-                      for (i = 0; i < options->horizon; i++)
-                        for (j = 0; j < nvars; j++)
-                          out_buf[(k) + ((i) + (j)*options->horizon)*options->num_percentiles] = ElementM(forecast, i, j);
-                    mxFree(forecast);
-                  }
-            }
-          mxFree(histogram);
-        }
-    }
-  catch (const char *s)
-    {
-      DYN_MEX_FUNC_ERR_MSG_TXT("Exception Thrown in Forecast: \n");
-    }
-
-  mxFree(p);
-  mxFree(model);
-  plhs[0] = mxCreateDoubleScalar(0);
-}
-
-#endif
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index 9702bf9b30f6f6df7a47a43d514194ee1974df85..56c5df7a7c04505fa88182ba2fa929b59e110df4 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -1093,6 +1093,13 @@ void
 MSSBVARForecastStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
 {
   mod_file_struct.bvar_present = true;
+
+  if (options_list.num_options.find("ms.regimes") != options_list.num_options.end())
+    if (options_list.num_options.find("ms.regime") != options_list.num_options.end())
+      {
+        cerr << "ERROR: You may only pass one of regime and regimes to ms_forecast" << endl;
+        exit(EXIT_FAILURE);
+      }
 }
 
 void
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 43781f5618572663dce744e9703e05d75dc624bb..abc9e1351ad8deacfdf8a5fe0d7aac461950815f 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -1652,13 +1652,15 @@ ms_forecast_option : o_output_file_tag
                    | o_file_tag
                    | o_simulation_file_tag
                    | o_data_obs_nbr
-                   | o_no_error_bands
                    | o_error_band_percentiles
                    | o_shock_draws
                    | o_shocks_per_parameter
                    | o_thinning_factor
                    | o_free_parameters
                    | o_median
+                   | o_regime
+                   | o_regimes
+                   | o_parameter_uncertainty
                    ;
 
 ms_forecast_options_list : ms_forecast_option COMMA ms_forecast_options_list