diff --git a/doc/dynare.texi b/doc/dynare.texi
index b3a25bd41305b90f3e6185f0046fc4428f7ef01c..e3049ec56b704df72e977600fd0eff2759a0267d 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -5858,16 +5858,13 @@ The forecast horizon. Default: @code{12}
 
 @item filtered_probabilities
 @anchor{filtered_probabilities} Uses filtered probabilities at the end
-of the sample as initial conditions for regime probabilities. Default:
-@code{off}
-
-@item no_error_bands
-@anchor{no_error_bands} Do not output error bands. Default: @code{off}
-(@i{i.e.} output error bands)
+of the sample as initial conditions for regime probabilities. Only one
+of @code{filtered_probabilities}, @code{regime} and @code{regimes} may
+be passed. Default: @code{off}
 
 @item error_band_percentiles = [@var{DOUBLE1} @dots{}]
 @anchor{error_band_percentiles} The percentiles to compute. Default:
-@code{[0.16 0.50 0.84]}. If @code{no_error_bands} is passed, the default
+@code{[0.16 0.50 0.84]}. If @code{median} is passed, the default
 is @code{[0.5]}
 
 @item shock_draws = @var{INTEGER}
@@ -5886,11 +5883,25 @@ draws in posterior draws file are used. Default: @code{1}
 @anchor{free_parameters} A vector of free parameters to initialize theta
 of the model. Default: use estimated parameters
 
-@item median
-@anchor{median}
+@item parameter_uncertainty
+@anchor{parameter_uncertainty} Calculate IRFs under parameter
+uncertainty. Requires that @command{ms_simulation} has been
+run. Default: @code{off}
 
-A shortcut to setting @code{error_band_percentiles=[0.5]}. Default:
-@code{off}
+@item regime = @var{INTEGER}
+@anchor{regime} Given the data and model parameters, what is the ergodic
+probability of being in the specified regime. Only one of
+@code{filtered_probabilities}, @code{regime} and @code{regimes} may be
+passed. Default: @code{off}
+
+@item regimes
+@anchor{regimes} Describes the evolution of regimes. Only one of
+@code{filtered_probabilities}, @code{regime} and @code{regimes} may be
+passed. Default: @code{off}
+
+@item median
+@anchor{median} A shortcut to setting
+@code{error_band_percentiles=[0.5]}. Default: @code{off}
 
 @end table
 
@@ -5923,7 +5934,8 @@ while data files are contained in @code{<output_file_tag/Forecast>}.
 The number of data points included in the output. Default: @code{0}
 
 @item no_error_bands
-@xref{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}.
diff --git a/matlab/ms-sbvar/initialize_ms_sbvar_options.m b/matlab/ms-sbvar/initialize_ms_sbvar_options.m
index 0afb1fc0bc0bb2ab86d309fe05f05b89351c42b7..3937641e8dce8dae9869758bc8844a242e8234d0 100644
--- a/matlab/ms-sbvar/initialize_ms_sbvar_options.m
+++ b/matlab/ms-sbvar/initialize_ms_sbvar_options.m
@@ -12,7 +12,7 @@ function options_=initialize_ms_sbvar_options(M_, options_)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2011 Dynare Team
+% Copyright (C) 2011-2012 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -100,12 +100,14 @@ options_.ms.real_time_smoothed_probabilities = 0;
 % irf
 options_.ms.horizon = 12;
 options_.ms.filtered_probabilities = 0;
-options_.ms.error_bands = 1;
 options_.ms.percentiles = [.16 .5 .84];
 options_.ms.parameter_uncertainty = 0;
 options_.ms.shock_draws = 10000;
 options_.ms.shocks_per_parameter = 10;
 options_.ms.median = 0;
+options_.ms.regime = 0;
+options_.ms.regimes = 0;
 % forecast
+options_.ms.error_bands = 1;
 options_.ms.forecast_data_obs = 0;
 end
diff --git a/matlab/ms-sbvar/ms_irf.m b/matlab/ms-sbvar/ms_irf.m
index c233395e9ae95a22f2ddd783f34456064d4f8f22..dcf89d3e3210cc8c271014c9367d908476a07a90 100644
--- a/matlab/ms-sbvar/ms_irf.m
+++ b/matlab/ms-sbvar/ms_irf.m
@@ -15,7 +15,7 @@ function [options_, oo_]=ms_irf(varlist, M_, options_, oo_)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2011 Dynare Team
+% Copyright (C) 2011-2012 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -35,50 +35,91 @@ function [options_, oo_]=ms_irf(varlist, M_, options_, oo_)
 disp('MS-SBVAR Impulse Response Function');
 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_, 'irf')
+clean_ms_irf_files(options_.ms.output_file_tag);
 irfdir = [options_.ms.output_file_tag filesep 'IRF'];
 create_dir(irfdir);
 
-opt = { ...
-    {'file_tag', options_.ms.file_tag}, ...
-    {'seed', options_.DynareRandomStreams.seed}, ...
-    {'horizon', options_.ms.horizon}, ...
-    {'filtered', options_.ms.filtered_probabilities}, ...
-    {'error_bands', options_.ms.error_bands}, ...
-    {'percentiles', options_.ms.percentiles}, ...
-    {'thin', options_.ms.thinning_factor}
-    };
+% setup command line options
+opt = ['-ir -seed ' num2str(options_.DynareRandomStreams.seed)];
+opt = [opt ' -do ' irfdir];
+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)];
 
-if options_.ms.median
-    opt = [opt(:)' {{'median'}}];
+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)];
+elseif options_.ms.filtered_probabilities
+    opt = [opt ' -filtered'];
 end
 
-[err, irf] = mex_ms_irf([opt(:)', {{'free_parameters', oo_.ms.maxparams}, {'shocks_per_parameter', options_.ms.shock_draws}}]);
-mexErrCheck('mex_ms_irf ergodic ', err);
-plot_ms_irf(M_,options_,irf,options_.varobs,'Ergodic Impulse Responses',varlist);
+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
 
-[err, regime_irfs] = mex_ms_irf([opt(:)', {{'free_parameters',oo_.ms.maxparams}, {'shocks_per_parameter', options_.ms.shock_draws}, {'regimes'}}]);
-mexErrCheck('mex_ms_irf ergodic regimes ',err);
-for i=1:size(regime_irfs,1)
-    plot_ms_irf(M_,options_,squeeze(regime_irfs(i,:,:,:)),options_.varobs,['Ergodic ' ...
-                        'Impulse Responses State ' int2str(i)],varlist);
+percentiles_size = 0;
+if options_.ms.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
-save([irfdir filesep 'ergodic_irf.mat'], 'irf', 'regime_irfs');
 
-if exist(options_.ms.mh_file,'file') > 0
-    [err, irf] = mex_ms_irf([opt(:)', {{'shocks_per_parameter', options_.ms.shocks_per_parameter}, ...
-        {'parameter_uncertainty'},{'simulation_file',options_.ms.mh_file}}]);
-    mexErrCheck('mex_ms_irf bayesian ',err);
-    plot_ms_irf(M_,options_,irf,options_.varobs,'Impulse Responses with Parameter Uncertainty',varlist);
-    
-    [err, regime_irfs] = mex_ms_irf([opt(:)', {{'shocks_per_parameter', options_.ms.shocks_per_parameter}, ...
-        {'simulation_file',options_.ms.mh_file},{'parameter_uncertainty'},{'regimes'}}]);
-    mexErrCheck('mex_ms_irf bayesian regimes ',err);
-    for i=1:size(regime_irfs,1)
-        plot_ms_irf(M_,options_,squeeze(regime_irfs(i,:,:,:)),options_.varobs,['Impulse ' ...
-                            'Responses with Parameter Uncertainty State ' int2str(i)],varlist);
+% irf
+[err] = ms_sbvar_command_line(opt);
+mexErrCheck('ms_irf',err);
+
+% Plot IRFs
+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
-    save([irfdir filesep 'bayesian_irf.mat'], 'irf', 'regime_irfs');
+
+    for state_i=1:n_states
+        irf_data = load([irfdir filesep 'ir_percentiles_regime_' ...
+            num2str(state_i-1) '_' options_.ms.output_file_tag ...
+            '.out'], '-ascii');
+        irf_data = reshape_ascii_irf_data(M_.endo_nbr, percentiles_size, ...
+            options_.ms.horizon, irf_data);
+        save([irfdir filesep 'irf_state_' num2str(state_i-1)], 'irf_data');
+        plot_ms_irf(M_,options_,irf_data,options_.varobs, ...
+            ['Impulse Responses, State ...' num2str(state_i)], varlist);
+    end
+else
+    if options_.ms.regime
+        irf_data = load([irfdir filesep 'ir_percentiles_regime_' ...
+            num2str(options_.ms.regime-1) '_' options_.ms.output_file_tag ...
+            '.out'], '-ascii');
+        irf_title = ['Impulse Response, State ' num2str(options_.ms.regime)];
+        save_filename = ['irf_regime_' num2str(options_.ms.regime-1)];
+    elseif options_.ms.filtered_probabilities
+        irf_data = load([irfdir filesep 'ir_percentiles_filtered_' ...
+            options_.ms.output_file_tag '.out'], '-ascii');
+        irf_title = 'Impulse Response Filtered';
+        save_filename = 'irf';
+    else
+        irf_data = load([irfdir filesep 'ir_percentiles_ergodic_' ...
+            options_.ms.output_file_tag '.out'], '-ascii');
+        irf_title = 'Impulse Response Ergodic';
+        save_filename = 'irf';
+    end
+
+    irf_data = reshape_ascii_irf_data(M_.endo_nbr, percentiles_size, ...
+        options_.ms.horizon, irf_data);
+    save([irfdir filesep save_filename], 'irf_data');
+    plot_ms_irf(M_, options_, irf_data, options_.varobs, irf_title, varlist);
 end
 end
diff --git a/matlab/ms-sbvar/reshape_ascii_irf_data.m b/matlab/ms-sbvar/reshape_ascii_irf_data.m
new file mode 100644
index 0000000000000000000000000000000000000000..3d9629b224cc3cdbd2f90fc54517c0065da5b72f
--- /dev/null
+++ b/matlab/ms-sbvar/reshape_ascii_irf_data.m
@@ -0,0 +1,44 @@
+function irf_data=reshape_ascii_irf_data(endo_nbr, psize, horizon, ascii_data)
+% function irf_data=reshape_ascii_irf_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
+%    irf_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
+    irf_data = ascii_data;
+    return
+end
+
+irf_data = zeros(psize, horizon, endo_nbr*endo_nbr);
+for i=1:endo_nbr*endo_nbr
+    for j=1:psize
+        irf_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 95cfb83ac30d50564989af7ed1c49d06cfcce4c1..1f703bbed7ff4e7ac022bd3163c20d3e85fb3e36 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_irf 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_forecast mex_ms_variance_decomposition
 
 DWSWITCHDIR = $(top_srcdir)/../../../ms-sbvar/switch_dw
 DWUTILITIESDIR = $(top_srcdir)/../../../ms-sbvar/utilities_dw
@@ -75,31 +75,7 @@ nodist_mex_ms_convert_free_parameters_SOURCES = \
 	$(MSMEXSRCDIR)/mex_ms_sbvar.cc \
   $(MSMEXSRCDIR)/modify_for_mex.cc \
 	$(DWUTILS)
-	
-	
-nodist_mex_ms_irf_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_irf.cc \
-	$(DWUTILS)
-	
+
 nodist_mex_ms_forecast_SOURCES = \
   $(DWSWITCHDIR)/switching/dw_switch.c \
   $(DWSWITCHDIR)/switching/dw_switchio.c \
diff --git a/mex/sources/ms-sbvar/mex_ms_irf.cc b/mex/sources/ms-sbvar/mex_ms_irf.cc
deleted file mode 100644
index ada45048710d0ac7ef223a13439a7e38323cc74c..0000000000000000000000000000000000000000
--- a/mex/sources/ms-sbvar/mex_ms_irf.cc
+++ /dev/null
@@ -1,190 +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"
-
-extern "C"  {
-#include "modify_for_mex.h"
-#include "switch.h"
-#include "switchio.h"
-#include "VARio.h"
-#include "dw_rand.h"
-#include "dw_histogram.h"
-#include "sbvar_impulse_responses.h"
-}
-
-void
-mexFunction(int nlhs, mxArray *plhs[],
-            int nrhs, const mxArray *prhs[])
-{
-  double *out_buf;
-  int i, j, k, s, nfree, nstates, nvars, npre;
-
-  TStateModel *model;
-  SbvarOption *options = (SbvarOption *) NULL;
-  T_VAR_Parameters *p = (T_VAR_Parameters *) NULL;
-  TMatrixHistogram *histogram = (TMatrixHistogram *) NULL;
-  TMatrix ir;
-
-  int type = F_FREE, ergodic = 1;
-
-  // Check the left hand and right hand side arguments to make sure they conform
-  if (nrhs != 1)
-    DYN_MEX_FUNC_ERR_MSG_TXT("ms_irf takes one cell array as an input argument.");
-  if (nlhs != 2)
-    DYN_MEX_FUNC_ERR_MSG_TXT("You must specify two output arguments.");
-
-  // copy the string data from prhs[0] into a C string input_ buf.
-  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);
-
-  // 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*nvar)
-        mwSize dims[4] = {nstates, options->num_percentiles, options->horizon, (nvars*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*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*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*nvars)};
-      plhs[1] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL);
-      out_buf = mxGetPr(plhs[1]);
-    }
-
-  // Use filter probabilities?
-  if (options->filtered_probabilities)
-    ergodic = 0;
-
-  try
-    {
-      if (options->regimes)
-        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++)
-                  if (ir = IRHistogramToPercentile((TMatrix) NULL, options->horizon, options->percentiles[k], histogram, model))
-                    {
-                      // 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*nvars); j++)
-                            out_buf[(s) + ((i) + (j)*options->horizon)*nstates] = ElementM(ir, i, j);
-                      else
-                        for (i = 0; i < options->horizon; i++)
-                          for (j = 0; j < (nvars*nvars); j++)
-                            out_buf[(s) + ((k) + ((i) + (j)*options->horizon)*options->num_percentiles)*nstates] = ElementM(ir, i, j);
-                      mxFree(ir);
-                    }
-              }
-            mxFree(histogram);
-          }
-
-      else if (options->regime >= 0)
-        {
-          s = options->regime;
-          printf("Constructing percentiles for impulse responses - regime %d\n", s);
-          if (histogram = impulse_response_percentile_regime(options->simulation_file, options->thin, options->regime, options->horizon, model, type))
-            {
-              for (k = 0; k < options->num_percentiles; k++)
-                if (ir = IRHistogramToPercentile((TMatrix) NULL, options->horizon, options->percentiles[k], histogram, model))
-                  {
-                    // 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*nvars); j++)
-                          out_buf[(i) + (j)*options->horizon] = ElementM(ir, i, j);
-                    else
-                      for (i = 0; i < options->horizon; i++)
-                        for (j = 0; j < (nvars*nvars); j++)
-                          out_buf[(k) + ((i) + (j)*options->horizon)*options->num_percentiles] = ElementM(ir, i, j);
-                    mxFree(ir);
-                  }
-            }
-          mxFree(histogram);
-        }
-      else
-        {
-          printf(ergodic ? "Constructing percentiles for ergodic impulse responses - %d draws of regimes per posterior value\n"
-                 : "Constructing percentiles for filtered impulse responses - %d draws of regimes per posterior value\n", options->shocks);
-          if (histogram = impulse_response_percentile_ergodic(options->shocks, options->simulation_file, options->thin, ergodic, options->horizon, model, type))
-            {
-              for (k = 0; k < options->num_percentiles; k++)
-                if (ir = IRHistogramToPercentile((TMatrix) NULL, options->horizon, options->percentiles[k], histogram, model))
-                  {
-                    // 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*nvars); j++)
-                          out_buf[(i) + (j)*options->horizon] = ElementM(ir, i, j);
-                    else
-                      for (i = 0; i < options->horizon; i++)
-                        for (j = 0; j < (nvars*nvars); j++)
-                          out_buf[(k) + ((i) + (j)*options->horizon)*options->num_percentiles] = ElementM(ir, i, j);
-                    mxFree(ir);
-                  }
-            }
-          mxFree(histogram);
-        }
-    }
-  catch (const char *s)
-    {
-      DYN_MEX_FUNC_ERR_MSG_TXT("Exception Thrown in IRF: \n");
-    }
-
-  mxFree(model);
-  mxFree(p);
-  plhs[0] = mxCreateDoubleScalar(0);
-}
-
-#endif
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index a8ac25d411714594268ab8e4ef2468a930320606..9702bf9b30f6f6df7a47a43d514194ee1974df85 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -1048,6 +1048,31 @@ void
 MSSBVARIrfStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
 {
   mod_file_struct.bvar_present = true;
+
+  bool regime_present = false;
+  bool regimes_present = false;
+  bool filtered_probabilities_present = false;
+
+  OptionsList::num_options_t::const_iterator it = options_list.num_options.find("ms.regimes");
+  if (it != options_list.num_options.end())
+    regimes_present = true;
+
+  it = options_list.num_options.find("ms.regime");
+  if (it != options_list.num_options.end())
+    regime_present = true;
+
+  it = options_list.num_options.find("ms.filtered_probabilities");
+  if (it != options_list.num_options.end())
+    filtered_probabilities_present = true;
+
+  if ((filtered_probabilities_present && regime_present) ||
+      (filtered_probabilities_present && regimes_present) ||
+      (regimes_present && regime_present))
+      {
+        cerr << "ERROR: You may only pass one of regime, regimes and "
+             << "filtered_probabilities to ms_irf" << endl;
+        exit(EXIT_FAILURE);
+      }
 }
 
 void
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 40857352d314936f48c81608ed7f7e783bc19946..43781f5618572663dce744e9703e05d75dc624bb 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2011 Dynare Team
+ * Copyright (C) 2003-2012 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -122,9 +122,9 @@ class ParsingDriver;
 %token PRINT PRIOR_MC PRIOR_TRUNC PRIOR_MODE PRIOR_MEAN POSTERIOR_MODE POSTERIOR_MEAN POSTERIOR_MEDIAN PRUNING
 %token <string_val> QUOTED_STRING
 %token QZ_CRITERIUM FULL DSGE_VAR DSGE_VARLAG DSGE_PRIOR_WEIGHT
-%token RELATIVE_IRF REPLIC RPLOT SAVE_PARAMS_AND_STEADY_STATE
+%token RELATIVE_IRF REPLIC RPLOT SAVE_PARAMS_AND_STEADY_STATE PARAMETER_UNCERTAINTY
 %token SHOCKS SHOCK_DECOMPOSITION SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED SMOOTHER STACK_SOLVE_ALGO STEADY_STATE_MODEL SOLVE_ALGO
-%token STDERR STEADY STOCH_SIMUL SYLVESTER
+%token STDERR STEADY STOCH_SIMUL SYLVESTER REGIMES REGIME
 %token TEX RAMSEY_POLICY PLANNER_DISCOUNT DISCRETIONARY_POLICY
 %token <string_val> TEX_NAME
 %token UNIFORM_PDF UNIT_ROOT_VARS USE_DLL USEAUTOCORR GSA_SAMPLE_FILE
@@ -1674,15 +1674,17 @@ ms_forecast : MS_FORECAST ';'
 ms_irf_option : o_output_file_tag
               | o_file_tag
               | o_simulation_file_tag
+              | o_parameter_uncertainty
               | o_horizon
               | o_filtered_probabilities
-              | 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
               ;
 
 ms_irf_options_list : ms_irf_option COMMA ms_irf_options_list
@@ -2316,6 +2318,7 @@ o_use_mean_center : USE_MEAN_CENTER { driver.option_num("ms.use_mean_center","1"
 o_proposal_type : PROPOSAL_TYPE EQUAL INT_NUMBER { driver.option_num("ms.proposal_type",$3); }
 o_proposal_lower_bound : PROPOSAL_LOWER_BOUND EQUAL signed_number { driver.option_num("ms.proposal_lower_bound",$3); }
 o_proposal_upper_bound : PROPOSAL_UPPER_BOUND EQUAL signed_number { driver.option_num("ms.proposal_upper_bound",$3); }
+o_parameter_uncertainty : PARAMETER_UNCERTAINTY { driver.option_num("ms.parameter_uncertainty","1"); };
 o_horizon : HORIZON EQUAL INT_NUMBER { driver.option_num("ms.horizon",$3); };
 o_filtered_probabilities : FILTERED_PROBABILITIES { driver.option_num("ms.filtered_probabilities","1"); };
 o_real_time_smoothed : REAL_TIME_SMOOTHED { driver.option_num("ms.real_time_smoothed_probabilities","1"); };
@@ -2325,6 +2328,8 @@ o_shock_draws : SHOCK_DRAWS EQUAL INT_NUMBER { driver.option_num("ms.shock_draws
 o_shocks_per_parameter : SHOCKS_PER_PARAMETER EQUAL INT_NUMBER { driver.option_num("ms.shocks_per_parameter",$3); };
 o_free_parameters : FREE_PARAMETERS EQUAL vec_value { driver.option_num("ms.free_parameters",$3); };
 o_median : MEDIAN { driver.option_num("ms.median","1"); };
+o_regimes : REGIMES { driver.option_num("ms.regimes","1"); };
+o_regime : REGIME EQUAL INT_NUMBER { driver.option_num("ms.regime",$3); };
 o_data_obs_nbr : DATA_OBS_NBR EQUAL INT_NUMBER { driver.option_num("ms.forecast_data_obs",$3); };
 
 range : symbol ':' symbol
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index a5010b5bdc198b1b411b8b09f1c68e324504268a..54c42442fde7f4f889b44535804eba50914dfdc5 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2011 Dynare Team
+ * Copyright (C) 2003-2012 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -375,6 +375,7 @@ string eofbuff;
 <DYNARE_STATEMENT>no_create_init {return token::NO_CREATE_INIT;};
 <DYNARE_STATEMENT>simulation_file_tag {return token::SIMULATION_FILE_TAG;};
 <DYNARE_STATEMENT>horizon {return token::HORIZON;}
+<DYNARE_STATEMENT>parameter_uncertainty {return token::PARAMETER_UNCERTAINTY;}
 <DYNARE_STATEMENT>no_error_bands {return token::NO_ERROR_BANDS;}
 <DYNARE_STATEMENT>error_band_percentiles {return token::ERROR_BAND_PERCENTILES;}
 <DYNARE_STATEMENT>shock_draws {return token::SHOCK_DRAWS;}
@@ -382,6 +383,8 @@ string eofbuff;
 <DYNARE_STATEMENT>thinning_factor {return token::THINNING_FACTOR;}
 <DYNARE_STATEMENT>free_parameters {return token::FREE_PARAMETERS;}
 <DYNARE_STATEMENT>median {return token::MEDIAN;}
+<DYNARE_STATEMENT>regime {return token::REGIME;}
+<DYNARE_STATEMENT>regimes {return token::REGIMES;}
 <DYNARE_STATEMENT>data_obs_nbr {return token::DATA_OBS_NBR;}
 <DYNARE_STATEMENT>filtered_probabilities {return token::FILTERED_PROBABILITIES;}
 <DYNARE_STATEMENT>real_time_smoothed {return token::REAL_TIME_SMOOTHED;}
diff --git a/tests/ms-sbvar/test_ms_variances.mod b/tests/ms-sbvar/test_ms_variances.mod
index 6f6d2477bc621cec93db8737ec54ecaad502c504..034624ffeedee0a4dad9bb1962c2d240257a3d79 100644
--- a/tests/ms-sbvar/test_ms_variances.mod
+++ b/tests/ms-sbvar/test_ms_variances.mod
@@ -23,6 +23,6 @@ ms_estimation(datafile=data
 ms_simulation(mh_replic=1000);
 ms_compute_mdd;
 ms_compute_probabilities;
-ms_irf;
+ms_irf(parameter_uncertainty);
 ms_forecast;
 ms_variance_decomposition;