diff --git a/doc/dynare.texi b/doc/dynare.texi
index 944c952568e254cea4f46287f77abb9bd165efa9..7f97ab7e40091c369e88d14c866d3859970af791 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -5854,7 +5854,7 @@ model. Output @code{.eps} files are contained in
 @xref{simulation_file_tag}.
 
 @item horizon = @var{INTEGER}
-The forecast horizon. Default: @code{12}
+@anchor{horizon} The forecast horizon. Default: @code{12}
 
 @item filtered_probabilities
 @anchor{filtered_probabilities} Uses filtered probabilities at the end
@@ -5990,12 +5990,15 @@ are contained in @code{<output_file_tag/Variance_Decomposition>}.
 @item simulation_file_tag = @var{FILENAME}
 @xref{simulation_file_tag}.
 
+@item horizon = @var{INTEGER}
+@xref{horizon}.
+
 @item filtered_probabilities
 @xref{filtered_probabilities}.
 
 @item no_error_bands
-Do not output error bands. Default: @code{off}
-(@i{i.e.} output error bands)
+Do not output percentile error bands (@i{i.e.} compute mean). Default:
+@code{off} (@i{i.e.} output error bands)
 
 @item error_band_percentiles = [@var{DOUBLE1} @dots{}]
 @xref{error_band_percentiles}.
@@ -6012,9 +6015,15 @@ Do not output error bands. Default: @code{off}
 @item free_parameters = @var{NUMERICAL_VECTOR}
 @xref{free_parameters}.
 
-@item median
+@item parameter_uncertainty
+@xref{parameter_uncertainty}.
 
-@xref{median}.
+@item regime = @var{INTEGER}
+@xref{regime}.
+
+@item regimes
+
+@xref{regimes}.
 
 @end table
 
diff --git a/matlab/ms-sbvar/ms_variance_decomposition.m b/matlab/ms-sbvar/ms_variance_decomposition.m
index 5469465b48b3fe558c7d5243a60ecfeb2c89db28..8b9822f130a8a06f23d416931939043bc10c4fc0 100644
--- a/matlab/ms-sbvar/ms_variance_decomposition.m
+++ b/matlab/ms-sbvar/ms_variance_decomposition.m
@@ -34,46 +34,100 @@ function [options_, oo_]=ms_variance_decomposition(M_, options_, oo_)
 disp('MS-SBVAR Variance Decomposition');
 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_, 'variance_decomposition')
+clean_ms_variance_decomposition_files(options_.ms.output_file_tag);
 vddir = [options_.ms.output_file_tag filesep 'Variance_Decomposition'];
 create_dir(vddir);
 
-% NOTICE THAT VARIANCE DECOMPOSITION DEFAULTS TO USING THE MEAN, NOT MEDIAN OR BANDED
+% setup command line options
+opt = ['-variance_decomposition -seed ' num2str(options_.DynareRandomStreams.seed)];
+opt = [opt ' -do ' vddir];
+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 = {
-    {'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}, ...
-    {'mean'} ...
-    };
+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
+
+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
 
-if options_.ms.median
-    opt = [opt(:)' {{'median'}}];
+percentiles_size = 1;
+outfile = [vddir filesep 'var_decomp_mean_'];
+if options_.ms.error_bands
+    % error_bands / percentiles used differently by
+    % Dan's variance decomposition code
+    % no_error_bands => mean is computed
+    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
+    outfile = [vddir filesep 'var_decomp_percentiles_'];
 end
 
-[err, vd] = mex_ms_variance_decomposition([opt(:)', {{'free_parameters',oo_.ms.maxparams}, ...
-    {'shocks_per_parameter', options_.ms.shock_draws}}]);
-mexErrCheck('mex_ms_variance_decomposition ergodic ', err);
-plot_ms_variance_decomposition(M_,options_,vd, 'Ergodic Variance Decomposition',options_.graph_save_formats,options_.TeX);
+% variance_decomposition
+[err] = ms_sbvar_command_line(opt);
+mexErrCheck('ms_variance_decomposition',err);
 
-[err, regime_vd] = mex_ms_variance_decomposition([opt(:)', {{'free_parameters',oo_.ms.maxparams}, ...
-    {'shocks_per_parameter', options_.ms.shock_draws}, {'regimes'}}]);
-mexErrCheck('mex_ms_variance_decomposition ergodic regimes', err);
-save([vddir filesep 'ergodic_vd.mat'], 'vd', 'regime_vd');
+if options_.ms.regime || options_.ms.regimes
+    outfile = [outfile 'regime_'];
+    if options_.ms.regime
+        outfile = [outfile num2str(options_.ms.regime-1) ...
+            '_' options_.ms.output_file_tag '.out'];
+    end
+elseif options_.ms.filtered_probabilities
+    outfile = [outfile 'filtered_' options_.ms.output_file_tag '.out'];
+else
+    outfile = [outfile 'ergodic_' options_.ms.output_file_tag '.out'];
+end
 
-if exist(options_.ms.mh_file,'file') > 0
-    [err, vd] = mex_ms_variance_decomposition([opt(:)', {{'simulation_file',options_.ms.mh_file}, ...
-        {'shocks_per_parameter', options_.ms.shocks_per_parameter}, {'parameter_uncertainty'}}]);
-    mexErrCheck('mex_ms_variance_decomposition bayesian ', err);
+% Create plots
+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
 
-    [err, regime_vd] = mex_ms_variance_decomposition([opt(:)', {{'simulation_file',options_.ms.mh_file}, ...
-        {'shocks_per_parameter', options_.ms.shocks_per_parameter}, {'parameter_uncertainty'}, {'regimes'}}]);
-    mexErrCheck('mex_ms_variance_decomposition bayesian regimes ', err);
-    save([vddir filesep 'bayesian_vd.mat'], 'vd', 'regime_vd');
+    for state_i=1:n_states
+        vd_data = load([outfile num2str(state_i-1) '_' ...
+            options_.ms.output_file_tag '.out'], '-ascii');
+        vd_data = reshape_ascii_variance_decomposition_data( ...
+            M_.endo_nbr, percentiles_size, options_.ms.horizon, vd_data);
+        save([vddir filesep 'variance_decomposition_state_' num2str(state_i-1)], 'vd_data');
+        plot_ms_variance_decomposition(M_, options_, vd_data, ...
+            ['Variance Decomposition, State ...' num2str(state_i)], ...
+            options_.graph_save_formats, options_.TeX);
+    end
+else
+    if options_.ms.regime
+        vd_title = ['Variance Decomposition, State ' num2str(options_.ms.regime)];
+        save_filename = ['variance_decomposition_regime_' num2str(options_.ms.regime-1)];
+    else
+        save_filename = 'variance_decomposition';
+        if options_.ms.filtered_probabilities
+            vd_title = 'Variance Decomposition Filtered';
+        else
+            vd_title = 'Variance Decomposition Ergodic';
+        end
+    end
+    vd_data = load(outfile, '-ascii');
+    vd_data = reshape_ascii_variance_decomposition_data( ...
+        M_.endo_nbr, percentiles_size, options_.ms.horizon, vd_data);
+    save([vddir filesep save_filename], 'vd_data');
+    plot_ms_variance_decomposition(M_, options_, vd_data, ...
+        vd_title, options_.graph_save_formats, options_.TeX);
 end
 end
diff --git a/matlab/ms-sbvar/plot_ms_variance_decomposition.m b/matlab/ms-sbvar/plot_ms_variance_decomposition.m
index ba30b487f7a8f9a1b2d46ce7b462b39e8514212b..74b5e470d330783a71c7a61c8b484ede5d42e07d 100644
--- a/matlab/ms-sbvar/plot_ms_variance_decomposition.m
+++ b/matlab/ms-sbvar/plot_ms_variance_decomposition.m
@@ -34,6 +34,11 @@ function plot_ms_variance_decomposition(M_, options_, vd, title_, graph_save_for
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
+if length(size(vd)) == 3
+    plot_ms_variance_decomposition_error_bands(M_, options_, vd, title_);
+    return;
+end
+
     nvars = M_.endo_nbr;
     endo_names = M_.endo_names;
 
diff --git a/matlab/ms-sbvar/plot_ms_variance_decomposition_error_bands.m b/matlab/ms-sbvar/plot_ms_variance_decomposition_error_bands.m
new file mode 100644
index 0000000000000000000000000000000000000000..4ddecbd7b43285dab10b001c5cff9103ecebc6d2
--- /dev/null
+++ b/matlab/ms-sbvar/plot_ms_variance_decomposition_error_bands.m
@@ -0,0 +1,100 @@
+function plot_ms_variance_decomposition_error_bands(M_, options_, vddata, title_)
+% function plot_ms_variance_decomposition_error_bands(M_, options_, vddata, title_)
+% plots the variance decomposition with percentiles
+%
+% INPUTS
+%    M_:          (struct)    model structure
+%    options_:    (struct)    options
+%    vddata:      (matrix)    variance_decomposition (percentile, options_.ms.horizon, nvar
+%                             x nvar)
+%    title:       (string)    title
+%
+
+% 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/>.
+
+nvars = M_.endo_nbr;
+endo_names = M_.endo_names;
+var_list = endo_names(1:M_.orig_endo_nbr,:);
+
+names = {};
+tex_names = {};
+m = 1;
+for i = 1:size(var_list)
+    tmp = strmatch(var_list(i,:), endo_names, 'exact');
+    if isempty(tmp)
+        error([var_list(i,:) ' isn''t an endogenous variable'])
+    end
+    tex_name = deblank(M_.endo_names_tex(tmp,:));
+    if ~isempty(tex_name)
+        names{m} = deblank(var_list(i,:));
+        tex_names{m} = tex_name;
+        m = m + 1;
+    end
+end
+
+for i=1:M_.exo_nbr
+    tex_name = deblank(M_.exo_names_tex(i,:));
+    if ~isempty(tex_name)
+        names{m} = deblank(M_.exo_names(i,:));
+        tex_names{m} = tex_name;
+        m = m + 1;
+    end
+end
+
+dims = size(vddata);
+if length(dims) ~= 3
+    error('The variance decomposition matrix passed to be plotted does not appear to be the correct size');
+end
+num_percentiles = dims(1);
+
+if size(endo_names, 1) ~= nvars
+    error('The names passed are not the same length as the number of variables')
+end
+
+for s=1:nvars
+    shock = zeros(options_.ms.horizon, nvars, num_percentiles);
+    for n=1:num_percentiles
+        for i=1:nvars
+            shock(:,i,n) = vddata(n, :, ((i-1) + ((s-1)*nvars)+1));
+        end
+    end
+    plot_banded_vddata_for_shock(shock, nvars, endo_names, ...
+        deblank(endo_names(s,:)), title_, ...
+        [options_.ms.output_file_tag filesep 'Output' filesep 'Variance_Decomposition'], ...
+        options_, names, tex_names);
+end
+end
+
+function [fig] = plot_banded_vddata_for_shock(vddata, nvars, endo_names, ...
+    shock_name, title_, dirname, options_, names, tex_names)
+fig = figure('Name', title_);
+npercentiles = size(vddata,3);
+for k=1:nvars
+    subplot(ceil(sqrt(nvars)), ceil(sqrt(nvars)),k);
+    for nn=1:npercentiles
+        plot(vddata(:,k,nn))
+        hold on
+    end
+    hold off
+    disp([endo_names(k,:) ' contribution to ' shock_name]);
+    title([endo_names(k,:) ' contribution to ' shock_name]);
+end
+dyn_save_graph(dirname, [title_ ' ' shock_name], ...
+    options_.graph_save_formats, options_.TeX, names, tex_names, ...
+    [title_ ' ' shock_name]);
+end
diff --git a/matlab/ms-sbvar/reshape_ascii_variance_decomposition_data.m b/matlab/ms-sbvar/reshape_ascii_variance_decomposition_data.m
new file mode 100644
index 0000000000000000000000000000000000000000..a96833fd56153cbc613b41d874349467c8b88361
--- /dev/null
+++ b/matlab/ms-sbvar/reshape_ascii_variance_decomposition_data.m
@@ -0,0 +1,44 @@
+function vd_data=reshape_ascii_variance_decomposition_data(endo_nbr, psize, horizon, ascii_data)
+% function vd_data=reshape_ascii_vd_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
+%    vd_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
+    vd_data = ascii_data;
+    return
+end
+
+vd_data = zeros(psize, horizon, endo_nbr*endo_nbr);
+for i=1:endo_nbr*endo_nbr
+    for j=1:psize
+        vd_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 6270789b224366148fd7a810f891fcee55c31b77..219bf5d127a0b1b57920dc9630450394b52b9d3b 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_variance_decomposition
+noinst_PROGRAMS = ms_sbvar_create_init_file ms_sbvar_command_line mex_ms_convert_free_parameters
 
 DWSWITCHDIR = $(top_srcdir)/../../../ms-sbvar/switch_dw
 DWUTILITIESDIR = $(top_srcdir)/../../../ms-sbvar/utilities_dw
@@ -75,26 +75,3 @@ nodist_mex_ms_convert_free_parameters_SOURCES = \
 	$(MSMEXSRCDIR)/mex_ms_sbvar.cc \
 	$(MSMEXSRCDIR)/modify_for_mex.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 \
-	$(MSMEXSRCDIR)/mex_ms_variance_decomposition.cc \
-	$(DWUTILS)
diff --git a/mex/sources/ms-sbvar/mex_ms_variance_decomposition.cc b/mex/sources/ms-sbvar/mex_ms_variance_decomposition.cc
deleted file mode 100644
index 2ef5f5d8b02a668594acab31faad075a8266636a..0000000000000000000000000000000000000000
--- a/mex/sources/ms-sbvar/mex_ms_variance_decomposition.cc
+++ /dev/null
@@ -1,243 +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 "switch.h"
-#include "switchio.h"
-#include "VARio.h"
-#include "dw_rand.h"
-#include "dw_histogram.h"
-#include "sbvar_variance_decomposition.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 tvd;
-
-  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_variance_decomposition 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);
-
-  /* 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 variance decomposition - regime %d\n", s);
-
-              if (options->mean)
-                {
-                  if (tvd = variance_decomposition_mean_regime(options->simulation_file, options->thin, s, options->horizon, model, type))
-                    {
-                      NormalizeVarianceDecomposition(tvd, p);
-                      ReorderVarianceDecomposition(tvd, p);
-                      for (i = 0; i < options->horizon; i++)
-                        for (j = 0; j < (nvars*nvars); j++)
-                          out_buf[(s) + ((i) + (j)*options->horizon)*nstates] = ElementM(tvd, i, j);
-                    }
-                  mxFree(tvd);
-                }
-              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);
-                      for (k = 0; k < options->num_percentiles; k++)
-                        {
-                          MatrixPercentile(tvd, options->percentiles[k], histogram);
-                          ReorderVarianceDecomposition(tvd, p);
-                          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(tvd, 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(tvd, i, j);
-                        }
-                      mxFree(tvd);
-                    }
-                  mxFree(histogram);
-                }
-            }
-
-        }
-      else if (options->regime >= 0)
-        {
-          s = options->regime;
-          printf("Constructing variance decomposition - regime %d\n", s);
-          if (options->mean)
-            {
-              if (tvd = variance_decomposition_mean_regime(options->simulation_file, options->thin, s, options->horizon, model, type))
-                {
-                  NormalizeVarianceDecomposition(tvd, p);
-                  ReorderVarianceDecomposition(tvd, p);
-                  for (i = 0; i < options->horizon; i++)
-                    for (j = 0; j < (nvars*nvars); j++)
-                      out_buf[(s) + ((i) + (j)*options->horizon)*nstates] = ElementM(tvd, i, j);
-                }
-              mxFree(tvd);
-
-            }
-          else
-            {
-              if (histogram = variance_decomposition_percentiles_regime(options->simulation_file, options->thin, s, options->horizon, model, type))
-                {
-                  tvd = CreateMatrix(options->horizon, nvars*nvars);
-                  for (k = 0; k < options->num_percentiles; k++)
-                    {
-                      MatrixPercentile(tvd, options->percentiles[k], histogram);
-                      ReorderVarianceDecomposition(tvd, p);
-                      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(tvd, 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(tvd, i, j);
-                    }
-                  mxFree(tvd);
-                }
-              mxFree(histogram);
-            }
-        }
-      else
-        {
-
-          printf("Constructing variance decomposition - %d draws of regimes per posterior value\n", options->shocks);
-          if (options->mean)
-            {
-              if (tvd = variance_decomposition_mean(options->shocks, options->simulation_file, options->thin, ergodic, options->horizon, model, type))
-                {
-                  NormalizeVarianceDecomposition(tvd, p);
-                  ReorderVarianceDecomposition(tvd, p);
-                  for (i = 0; i < options->horizon; i++)
-                    for (j = 0; j < (nvars*nvars); j++)
-                      out_buf[(i) + (j)*options->horizon] = ElementM(tvd, i, j);
-                }
-              mxFree(tvd);
-            }
-          else
-            {
-
-              if (histogram = variance_decomposition_percentiles(options->shocks, options->simulation_file, options->thin, ergodic, options->horizon, model, type))
-                {
-
-                  tvd = CreateMatrix(options->horizon, nvars*nvars);
-                  for (k = 0; k < options->num_percentiles; k++)
-                    {
-                      MatrixPercentile(tvd, options->percentiles[k], histogram);
-                      ReorderVarianceDecomposition(tvd, p);
-                      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(tvd, 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(tvd, i, j);
-                    }
-                  mxFree(tvd);
-                }
-              mxFree(histogram);
-            }
-        }
-    }
-  catch (const char *s)
-    {
-      DYN_MEX_FUNC_ERR_MSG_TXT("Exception Thrown in Variance Decomposition: \n");
-    }
-
-  mxFree(p);
-  mxFree(model);
-  plhs[0] = mxCreateDoubleScalar(0);
-}
-
-#endif
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index 56c5df7a7c04505fa88182ba2fa929b59e110df4..525d8070f1e9e72a2669989caa5782cdf0751d2e 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -1119,6 +1119,31 @@ void
 MSSBVARVarianceDecompositionStatement::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_variance_decomposition" << endl;
+        exit(EXIT_FAILURE);
+      }
 }
 
 void
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index abc9e1351ad8deacfdf8a5fe0d7aac461950815f..f6b4ebf88da633efb453534cff435871a778be8e 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -1635,7 +1635,10 @@ ms_variance_decomposition_option : o_output_file_tag
                                  | o_shocks_per_parameter
                                  | o_thinning_factor
                                  | o_free_parameters
-                                 | o_median
+                                 | o_regime
+                                 | o_regimes
+                                 | o_parameter_uncertainty
+                                 | o_horizon
                                  ;
 
 ms_variance_decomposition_options_list : ms_variance_decomposition_option COMMA ms_variance_decomposition_options_list
diff --git a/tests/ms-sbvar/test_ms_variances.mod b/tests/ms-sbvar/test_ms_variances.mod
index 034624ffeedee0a4dad9bb1962c2d240257a3d79..d9bca5a987d4787f53cbe531b143b26c0f75e875 100644
--- a/tests/ms-sbvar/test_ms_variances.mod
+++ b/tests/ms-sbvar/test_ms_variances.mod
@@ -25,4 +25,4 @@ ms_compute_mdd;
 ms_compute_probabilities;
 ms_irf(parameter_uncertainty);
 ms_forecast;
-ms_variance_decomposition;
+ms_variance_decomposition(no_error_bands);