diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 1cfb81abd9655f35836cfc555db206dab48d525b..273a22cf541bb5f3be56bb311c4698a384ca0915 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -255,56 +255,12 @@ options_.initial_date.subperiod = 0;
 % discretionary policy
 options_.discretionary_policy = 0;
 
-% SWZ SBVAR
-% prepare ms sbvar
-options_.ms.coefficients_prior_hyperparameters = [1.0 1.0 0.1 1.2 1.0 1.0];
-options_.ms.freq = 1;
-options_.ms.initial_subperiod = 1;
-options_.ms.final_subperiod=4;
-options_.ms.nlags = 1;
-options_.ms.cross_restrictions = 0;
-options_.ms.contemp_reduced_form = 0;
-options_.ms.bayesian_prior = 1;
-options_.ms.alpha = 1.0;
-options_.ms.beta = 1.0;
-options_.ms.gsig2_lmd = 50^2;
-options_.ms.gsig2_lmdm = 50^2;
-options_.ms.lower_cholesky = 0;
-options_.ms.upper_cholesky = 0;
-% all mex functions
-options_.ms.output_file_tag = M_.fname;
-%simulate
-options_.ms.mh_replic = 10000; % default differs from Dan's code
-% mdd
-options_.ms.mdd_proposal_type = [3 0.1 0.9];
-% irf
-options_.ms.bayesian_irf = 0;
-options_.ms.thinning_factor = 1;
-options_.ms.shock_draws = 10000;
-options_.ms.shocks_per_parameter = 10;
-options_.ms.percentiles = [.16 .5 .84];
-options_.ms.mode_compute = 1;
-options_.ms.mode_file = 0;
-options_.ms.load_mh_file = 0;
-options_.ms.filtered_probabilities = 0;
-options_.ms.real_time_smoothed_probabilities = 0;
-options_.ms.irf_shocks_per_parameter = options_.ms.shocks_per_parameter;
-options_.ms.irf_shock_draws = options_.ms.shock_draws;
-options_.ms.irf_thinning_factor = 1;
-options_.ms.irf_filtered = 0;
-options_.ms.forecast_shocks_per_parameter = options_.ms.shocks_per_parameter;
-options_.ms.forecast_shock_draws = options_.ms.shock_draws;
-options_.ms.forecast_thinning_factor = 1;
-options_.ms.forecast_data_obs = 0;
-options_.ms.vd_shocks_per_parameter = options_.ms.shocks_per_parameter;
-options_.ms.vd_shock_draws = options_.ms.shock_draws;
-options_.ms.vd_thinning_factor = 1;
-options_.ms.vd_filtered = 0;
-
-
 % Shock decomposition
 options_.parameter_set = [];
 
+% MS Sbvar options
+options_ = initialize_ms_sbvar_options(M_, options_);
+
 % initialize persistent variables in priordens()
 priordens([],[],[],[],[],[],1);
 
diff --git a/matlab/initialize_ms_sbvar_options.m b/matlab/initialize_ms_sbvar_options.m
new file mode 100644
index 0000000000000000000000000000000000000000..6f33c075fba77638f711affedbf7460f519e0fdf
--- /dev/null
+++ b/matlab/initialize_ms_sbvar_options.m
@@ -0,0 +1,109 @@
+function options_=initialize_ms_sbvar_options(M_, options_)
+%function initialize_ms_sbvar_options()
+% sets ms sbvar options back to their default values
+%
+% INPUTS
+%    M_
+%    options_
+%
+% OUTPUTS
+%    options_
+%
+% SPECIAL REQUIREMENTS
+%    none
+
+% 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/>.
+
+% MS SBVAR
+% all mex functions
+options_.ms.output_file_tag = M_.fname;
+% prepare ms sbvar & estimation
+%options_.ms.mode_file
+options_.ms.coefficients_prior_hyperparameters = [1.0 1.0 0.1 1.2 1.0 1.0];
+options_.ms.freq = 1;
+options_.ms.initial_subperiod = 1;
+options_.ms.final_subperiod=4;
+%options_.ms_varlist
+options_.ms.nlags = 1;
+options_.ms.cross_restrictions = 0;
+options_.ms.contemp_reduced_form = 0;
+options_.ms.bayesian_prior = 1;
+options_.ms.alpha = 1.0;
+options_.ms.beta = 1.0;
+options_.ms.gsig2_lmd = 50^2;
+options_.ms.gsig2_lmdm = 50^2;
+options_.ms.lower_cholesky = 0;
+options_.ms.upper_cholesky = 0;
+if isfield(options_.ms,'initial_year')
+    options_.ms = rmfield(options_.ms,'initial_year');
+end
+if isfield(options_.ms,'final_year')
+    options_.ms = rmfield(options_.ms,'final_year');
+end
+if isfield(options_.ms,'datafile')
+    options_.ms = rmfield(options_.ms,'datafile');
+end
+if isfield(options_.ms,'varlist')
+    options_.ms = rmfield(options_.ms,'varlist');
+end
+% estimation
+options_.ms.convergence_starting_value = 1e-3;
+options_.ms.convergence_ending_value = 1e-6;
+options_.ms.convergence_increment_value = 0.1;
+options_.ms.max_iterations_starting_value = 50;
+options_.ms.max_iterations_increment_value = 2;
+options_.ms.max_block_iterations = 100;
+options_.ms.max_repeated_optimization_runs = 10;
+options_.ms.function_convergence_criterion = 0.1;
+options_.ms.parameter_convergence_criterion = 0.1;
+options_.ms.number_of_large_perturbations = 5;
+options_.ms.number_of_small_perturbations = 5;
+options_.ms.number_of_posterior_draws_after_perturbation = 1;
+options_.ms.max_number_of_stages = 20;
+options_.ms.random_function_convergence_criterion = 0.1;
+options_.ms.random_parameter_convergence_criterion = 0.1;
+% simulation
+options_.ms.mh_replic = 10000; % default differs from Dan's code
+options_.ms.drop = 0.1*options_.ms.mh_replic;
+options_.ms.thinning_factor = 1;
+options_.ms.adaptive_mh_draws = 30000;
+% mdd
+options_.ms.mdd_proposal_draws = 100000;
+options_.ms.mdd_use_mean_center = 0;
+options_.ms.mdd_proposal_type = [3 0.1 0.9];
+% probabilities
+options_.ms.filtered_probabilities = 0;
+options_.ms.real_time_smoothed_probabilities = 0;
+% irf
+options_.ms.horizon = 12;
+options_.ms.filtered_probabilities = 0;
+options_.ms.error_bands = 1;
+options_.ms.error_band_percentiles = [.16 .5 .84];
+options_.ms.parameter_uncertainty = 0;
+options_.ms.shock_draws = 10000;
+options_.ms.shocks_per_parameter = 10;
+options_.ms.median = 0;
+% forecast
+options_.ms.forecast_data_obs = 0;
+if isfield(options_.ms,'free_parameters')
+    options_.ms = rmfield(options_.ms,'free_parameters');
+end
+if isfield(options_.ms,'simulation_file')
+    options_.ms = rmfield(options_.ms,'simulation_file');
+end
+end
diff --git a/matlab/ms-sbvar/clean_ms_files.m b/matlab/ms-sbvar/clean_ms_files.m
index a19fee68e81d60523879e390d9ed50f8528c9e2e..dfba902a6b1f361e92696533067673a862f7b47a 100644
--- a/matlab/ms-sbvar/clean_ms_files.m
+++ b/matlab/ms-sbvar/clean_ms_files.m
@@ -35,14 +35,24 @@ delete_if_exist(['est_flat_header_' output_file_tag '.out']);
 delete_if_exist(['est_flat_' output_file_tag '.out']);
 delete_if_exist(['est_free_' output_file_tag '.out']);
 delete_if_exist(['est_intermediate_' output_file_tag '.out']);
+delete_if_exist(['filtered_' output_file_tag '.out']);
 delete_if_exist('g1.mat');
 delete_if_exist('H.dat');
 delete_if_exist(['init_' output_file_tag '.dat']);
 delete_if_exist(['matlab_' output_file_tag '.prn']);
+delete_if_exist(['mdd_t1_' output_file_tag '.out']);
+delete_if_exist(['proposal_t1_' output_file_tag '.out']);
+delete_if_exist(['mdd_t2_' output_file_tag '.out']);
+delete_if_exist(['proposal_t2_' output_file_tag '.out']);
 delete_if_exist(['mdd_t3_' output_file_tag '.out']);
 delete_if_exist(['proposal_t3_' output_file_tag '.out']);
+delete_if_exist(['mdd_t4_' output_file_tag '.out']);
+delete_if_exist(['proposal_t4_' output_file_tag '.out']);
+delete_if_exist(['mdd_t5_' output_file_tag '.out']);
+delete_if_exist(['proposal_t5_' output_file_tag '.out']);
 delete_if_exist(['simulation_info_' output_file_tag '.out']);
 delete_if_exist(['simulation_' output_file_tag '.out']);
+delete_if_exist(['smoothed_' output_file_tag '.out']);
 delete_if_exist([output_file_tag '_markov_file.dat']);
 end
 
diff --git a/matlab/ms-sbvar/ms_compute_mdd.m b/matlab/ms-sbvar/ms_compute_mdd.m
index 387c89737103c3dff9608ae5721e1d8ecb4b4e22..f77bbc6100cb743ec6a457694dd8732bdb7c62ab 100644
--- a/matlab/ms-sbvar/ms_compute_mdd.m
+++ b/matlab/ms-sbvar/ms_compute_mdd.m
@@ -1,12 +1,15 @@
-function oo_=ms_compute_mdd(options_,oo_)
+function [options_, oo_]=ms_compute_mdd(M_, options_, oo_)
 %function ms_compute_mdd()
-% calls ms sbvar mdd mex file
+% Compute marginal data density
 %
 % INPUTS
+%    M_
 %    options_
+%    oo_
 %
 % OUTPUTS
-%    none
+%    options_
+%    oo_
 %
 % SPECIAL REQUIREMENTS
 %    none
@@ -28,8 +31,28 @@ function oo_=ms_compute_mdd(options_,oo_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-compute_mdd = create_mdd_commandline(options_);
-[err] = ms_sbvar_command_line(compute_mdd);
+disp('Marginal Data Density');
+options_ = set_ms_estimation_flags_for_other_mex(options_);
+options_ = set_ms_simulation_flags_for_other_mex(options_);
+oo_ = set_oo_w_estimation_output(options_, oo_);
+
+% setup command line options
+opt = ['-mdd -seed ' num2str(options_.DynareRandomStreams.seed)];
+opt = [opt ' -ft ' options_.ms.output_file_tag];
+opt = [opt ' -d ' num2str(options_.ms.mdd_proposal_draws)];
+if options_.ms.mdd_use_mean_center
+    opt = [opt ' -use_mean'];
+end
+if size(options_.ms.mdd_proposal_type, 2) ~= 3
+    error('ERROR: mdd_proposal_type takes exactly 3 arguments');
+else
+    opt = [opt ' -pt ' num2str(options_.ms.mdd_proposal_type(1))];
+    opt = [opt ' -l '  num2str(options_.ms.mdd_proposal_type(2))];
+    opt = [opt ' -h '  num2str(options_.ms.mdd_proposal_type(3))];
+end
+
+% compute mdd
+[err] = ms_sbvar_command_line(opt);
 mexErrCheck('ms_sbvar_command_line mdd',err);
 
 % grab the muller/bridge log mdd from the output file
@@ -54,25 +77,5 @@ if exist(mdd_filename,'file')
     oo_.ms.mueller_log_mdd = muller_mdd;
     oo_.ms.bridged_log_mdd = bridge_mdd;
 end
-end
-
-function opt=create_mdd_commandline(options_)
-
-opt = '-mdd';
-opt = [opt set_flag_if_exists(options_.ms, 'output_file_tag')];
-opt = [opt set_flag_if_exists(options_.ms, 'load_mh_file ')];
-opt = [opt set_flag_if_exists(options_.ms, 'mdd_proposal_draws ')];
-opt = [opt set_flag_if_exists(options_.ms, 'mdd_use_mean_center')];
-
-if size(options_.ms.mdd_proposal_type, 2) ~= 3
-    error('ERROR: mdd_proposal_type takes exactly 3 arguments');
-else
-    opt = [opt ' -pt ' num2str(options_.ms.mdd_proposal_type(1))];
-    opt = [opt ' -l '  num2str(options_.ms.mdd_proposal_type(2))];
-    opt = [opt ' -h '  num2str(options_.ms.mdd_proposal_type(3))];
-end
-
-if options_.DynareRandomStreams.seed
-    opt = [' -seed' num2str(options_.DynareRandomStreams.seed)];
-end
+options_ = initialize_ms_sbvar_options(M_, options_);
 end
diff --git a/matlab/ms-sbvar/ms_compute_probabilities.m b/matlab/ms-sbvar/ms_compute_probabilities.m
index 42b6b07d8c88b72a73c1d49982a571e3c5473429..73c731b5612d80b469cce4d77d1a4f4d54f1c179 100644
--- a/matlab/ms-sbvar/ms_compute_probabilities.m
+++ b/matlab/ms-sbvar/ms_compute_probabilities.m
@@ -1,12 +1,15 @@
-function ms_compute_probabilities(options_)
+function [options_, oo_]=ms_compute_probabilities(M_, options_, oo_)
 %function ms_simulation()
-% calls ms sbvar probabilities mex file
+% Compute posterior mode regime probabilities
 %
 % INPUTS
+%    M_
 %    options_
+%    oo_
 %
 % OUTPUTS
-%    none
+%    options_
+%    oo_
 %
 % SPECIAL REQUIREMENTS
 %    none
@@ -28,20 +31,34 @@ function ms_compute_probabilities(options_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-compute_probabilities = create_probabilities_commandline(options_);
-[err] = ms_sbvar_command_line(compute_probabilities);
-mexErrCheck('ms_sbvar_command_line probabilities',err);
+disp('Compute Marginal Data Density');
+options_ = set_ms_estimation_flags_for_other_mex(options_);
+options_ = set_ms_simulation_flags_for_other_mex(options_);
+oo_ = set_oo_w_estimation_output(options_, oo_);
 
-end
+% setup command line options
+opt = ['-probabilities -seed ' num2str(options_.DynareRandomStreams.seed)];
+opt = [opt ' -ft ' options_.ms.output_file_tag];
 
-function opt=create_probabilities_commandline(options_)
+if options_.ms.filtered_probabilities
+    opt = [opt ' -filtered' ];
+    prob_out_file = ['filtered_' options_.ms.output_file_tag '.out'];
+elseif options_.ms.real_time_smoothed_probabilities
+    opt = [opt ' -real_time_smoothed' ];
+    prob_out_file = 0;
+else
+    opt = [opt ' -smoothed' ];
+    prob_out_file = ['smoothed_' options_.ms.output_file_tag '.out'];
+end
 
-opt = '-simulate';
-opt = [opt set_flag_if_exists(options_.ms, 'output_file_tag')];
-opt = [opt set_flag_if_exists(options_.ms, 'filtered_probabilities')];
-opt = [opt set_flag_if_exists(options_.ms, 'output_file_tag')];
+% compute probabilities
+[err] = ms_sbvar_command_line(opt);
+mexErrCheck('ms_sbvar_command_line probabilities',err);
 
-if options_.DynareRandomStreams.seed
-    opt = [' -seed' num2str(options_.DynareRandomStreams.seed)];
+% now we want to plot the probabilities for each chain
+if ischar(prob_out_file)
+    computed_probabilities = load(prob_out_file);
+    plot_ms_probabilities(options_, computed_probabilities);
 end
+options_ = initialize_ms_sbvar_options(M_, options_);
 end
diff --git a/matlab/ms-sbvar/ms_estimation.m b/matlab/ms-sbvar/ms_estimation.m
index 400e62773aebe4175e0fcf5787c5cce2cfffa365..2e1bedf0c983249fa092bcd924ecc0f99ff793f4 100644
--- a/matlab/ms-sbvar/ms_estimation.m
+++ b/matlab/ms-sbvar/ms_estimation.m
@@ -1,12 +1,15 @@
-function ms_estimation(options_)
+function [options_, oo_]=ms_estimation(M_, options_, oo_)
 %function ms_estimation()
-% calls ms sbvar estimation mex file
+% MS Sbvar Estimation
 %
 % INPUTS
+%    M_
 %    options_
+%    oo_
 %
 % OUTPUTS
-%    none
+%    options_
+%    oo_
 %
 % SPECIAL REQUIREMENTS
 %    none
@@ -28,45 +31,38 @@ function ms_estimation(options_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
+disp('Estimation');
+
 % cleanup old files
 clean_ms_files(options_.ms.output_file_tag);
 
 % general setup
 ms_sbvar_setup(options_);
 
-% estimation
-perform_estimation = create_estimation_commandline(options_);
-[err] = ms_sbvar_command_line(perform_estimation);
-mexErrCheck('ms_sbvar_command_line estimation',err);
-
-end
+% setup command line options
+opt = ['-estimate -seed ' num2str(options_.DynareRandomStreams.seed)];
+opt = [opt ' -ft ' options_.ms.output_file_tag];
+opt = [opt ' -cb ' num2str(options_.ms.convergence_starting_value)];
+opt = [opt ' -ce ' num2str(options_.ms.convergence_ending_value)];
+opt = [opt ' -ci ' num2str(options_.ms.convergence_increment_value)];
+opt = [opt ' -ib ' num2str(options_.ms.max_iterations_starting_value)];
+opt = [opt ' -ii ' num2str(options_.ms.max_iterations_increment_value)];
+opt = [opt ' -mb ' num2str(options_.ms.max_block_iterations)];
+opt = [opt ' -repeat_max ' num2str(options_.ms.max_repeated_optimization_runs)];
+opt = [opt ' -repeat_tol_obj ' num2str(options_.ms.function_convergence_criterion)];
+opt = [opt ' -repeat_tol_parms ' num2str(options_.ms.parameter_convergence_criterion)];
+opt = [opt ' -random ' num2str(options_.ms.number_of_large_perturbations)];
+opt = [opt ' -random_small ' num2str(options_.ms.number_of_small_perturbations)];
+opt = [opt ' -random_small_ndraws ' num2str(options_.ms.number_of_posterior_draws_after_perturbation)];
+opt = [opt ' -random_max ' num2str(options_.ms.max_number_of_stages)];
+opt = [opt ' -random_tol_obj ' num2str(options_.ms.random_function_convergence_criterion)];
+opt = [opt ' -random_tol_parms ' num2str(options_.ms.random_parameter_convergence_criterion)];
 
-function opt=create_estimation_commandline(options_)
-
-opt = '-estimate';
-opt = [opt set_flag_if_exists(options_.ms, 'output_file_tag')];
-opt = [opt set_flag_if_exists(options_.ms, 'convergence_starting_value')];
-opt = [opt set_flag_if_exists(options_.ms, 'convergence_ending_value')];
-opt = [opt set_flag_if_exists(options_.ms, 'convergence_increment_value')];
-opt = [opt set_flag_if_exists(options_.ms, 'max_iterations_starting_value')];
-opt = [opt set_flag_if_exists(options_.ms, 'max_iterations_increment_value')];
-opt = [opt set_flag_if_exists(options_.ms, 'max_block_iterations')];
-opt = [opt set_flag_if_exists(options_.ms, 'max_repeated_optimization_runs')];
-opt = [opt set_flag_if_exists(options_.ms, 'function_convergence_criterion')];
-opt = [opt set_flag_if_exists(options_.ms, 'parameter_convergence_criterion')];
-opt = [opt set_flag_if_exists(options_.ms, 'number_of_large_perturbations')];
-opt = [opt set_flag_if_exists(options_.ms, 'number_of_small_perturbations')];
-opt = [opt set_flag_if_exists(options_.ms, 'number_of_posterior_draws_after_perturbation')];
-opt = [opt set_flag_if_exists(options_.ms, 'max_number_of_stages')];
-opt = [opt set_flag_if_exists(options_.ms, 'random_function_convergence_criterion')];
-opt = [opt set_flag_if_exists(options_.ms, 'random_parameter_convergence_criterion')];
-
-if options_.DynareRandomStreams.seed
-    opt = [' -seed' num2str(options_.DynareRandomStreams.seed)];
-end
-
-% if ischar(options_.ms.mode_file) > 0 && exist(options_.ms.mode_file) > 0
-%     opt = [opt ' -fp ' options_.ms.mode_file];
-% end
+% estimation
+[err] = ms_sbvar_command_line(opt);
+mexErrCheck('ms_sbvar_command_line estimation', err);
 
+options_ = set_ms_estimation_flags_for_other_mex(options_);
+oo_ = set_oo_w_estimation_output(options_, oo_);
+options_ = initialize_ms_sbvar_options(M_, options_);
 end
diff --git a/matlab/ms-sbvar/ms_forecast.m b/matlab/ms-sbvar/ms_forecast.m
new file mode 100644
index 0000000000000000000000000000000000000000..f9d59546bcbf6cafadcbf4f1e2f4f93c0b951ea4
--- /dev/null
+++ b/matlab/ms-sbvar/ms_forecast.m
@@ -0,0 +1,67 @@
+function [options_, oo_]=ms_forecast(M_, options_, oo_)
+%function ms_forecast()
+% calls ms irf mex function
+%
+% INPUTS
+%    M_
+%    options_
+%    oo_
+%
+% OUTPUTS
+%    options_
+%    oo_
+%
+% SPECIAL REQUIREMENTS
+%    none
+
+% 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/>.
+
+disp('Forecasts');
+options_ = set_ms_estimation_flags_for_other_mex(options_);
+options_ = set_ms_simulation_flags_for_other_mex(options_);
+oo_ = set_oo_w_estimation_output(options_, oo_);
+
+opt = {options_.ms.output_file_tag, ...
+    'seed', options_.DynareRandomStreams.seed, ...
+    'horizon', options_.ms.horizon, ...
+    'data',options_.ms.forecast_data_obs ...
+    'error_bands', options_.ms.error_bands, ...
+    'percentiles', options_.ms.error_band_percentiles, ...
+    'thin', options_.ms.thinning_factor };
+
+[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_,forecast,'Forecast');
+
+[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([M_.fname '/' options_.ms.output_file_tag '_ergodic_forecast.mat'], 'forecast', 'regime_forecast');
+
+if exist(options_.ms.load_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.load_mh_file,'parameter_uncertainty');
+    mexErrCheck('mex_ms_forecast bayesian ', err);
+    plot_ms_forecast(M_,forecast,'Forecast w/ Parameter Uncertainty');
+
+    [err, regime_forecast] = mex_ms_forecast(opt{:},'free_parameters',oo_.ms.maxparams,'shocks_per_parameter', options_.ms.shocks_per_parameter, ...
+        'simulation_file',options_.ms.load_mh_file,'parameter_uncertainty','regimes');
+    mexErrCheck('mex_ms_forecast bayesian regimes ', err);
+    save([M_.fname '/' options_.ms.output_file_tag '_bayesian_forecast.mat'], 'forecast', 'regime_forecast');
+end
+options_ = initialize_ms_sbvar_options(M_, options_);
+end
diff --git a/matlab/ms-sbvar/ms_irf.m b/matlab/ms-sbvar/ms_irf.m
new file mode 100644
index 0000000000000000000000000000000000000000..42ebe0950ffe52f4f52fc8fcf8a43ee3990ee11b
--- /dev/null
+++ b/matlab/ms-sbvar/ms_irf.m
@@ -0,0 +1,67 @@
+function [options_, oo_]=ms_irf(M_, options_, oo_)
+%function ms_irf()
+% calls ms irf mex function
+%
+% INPUTS
+%    M_
+%    options_
+%    oo_
+%
+% OUTPUTS
+%    options_
+%    oo_
+%
+% SPECIAL REQUIREMENTS
+%    none
+
+% 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/>.
+
+disp('IRFs');
+options_ = set_ms_estimation_flags_for_other_mex(options_);
+options_ = set_ms_simulation_flags_for_other_mex(options_);
+oo_ = set_oo_w_estimation_output(options_, oo_);
+
+opt = {options_.ms.output_file_tag, ...
+    'seed', options_.DynareRandomStreams.seed, ...
+    'horizon', options_.ms.horizon, ...
+    'filtered', options_.ms.filtered_probabilities, ...
+    'error_bands', options_.ms.error_bands, ...
+    'percentiles', options_.ms.error_band_percentiles, ...
+    'thin', options_.ms.thinning_factor };
+
+[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_,irf,options_.varobs,'Ergodic Impulse Responses');
+
+[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);
+save([M_.fname '/' options_.ms.output_file_tag '_ergodic_irf.mat'], 'irf', 'regime_irfs');
+
+if exist(options_.ms.load_mh_file,'file') > 0
+    [err, irf] = mex_ms_irf(opt{:}, 'shocks_per_parameter', options_.ms.shocks_per_parameter, ...
+        'parameter_uncertainty','simulation_file',options_.ms.load_mh_file);
+    mexErrCheck('mex_ms_irf bayesian ',err);
+    plot_ms_irf(M_,irf,options_.varobs,'Impulse Responses w/ Parameter Uncertainty');
+    
+    [err, regime_irfs] = mex_ms_irf(opt{:}, 'shocks_per_parameter', options_.ms.shocks_per_parameter, ...
+        'simulation_file',options_.ms.load_mh_file,'parameter_uncertainty','regimes');
+    mexErrCheck('mex_ms_irf bayesian regimes ',err);
+    save([M_.fname '/' options_.ms.output_file_tag '_bayesian_irf.mat'], 'irf', 'regime_irfs');
+end
+options_ = initialize_ms_sbvar_options(M_, options_);
+end
diff --git a/matlab/ms-sbvar/ms_simulation.m b/matlab/ms-sbvar/ms_simulation.m
index adbc0cb5812cee0187cf581a3b70980647519b49..26e742346fc02cbe3bd82205429f300462ca6bcd 100644
--- a/matlab/ms-sbvar/ms_simulation.m
+++ b/matlab/ms-sbvar/ms_simulation.m
@@ -1,12 +1,15 @@
-function ms_simulation(options_)
+function [options_, oo_]=ms_simulation(M_, options_, oo_)
 %function ms_simulation()
-% calls ms sbvar simulation mex file
+% MS Sbvar Simulation
 %
 % INPUTS
+%    M_
 %    options_
+%    oo_
 %
 % OUTPUTS
-%    none
+%    options_
+%    oo_
 %
 % SPECIAL REQUIREMENTS
 %    none
@@ -28,22 +31,22 @@ function ms_simulation(options_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-perform_simulation = create_simulation_commandline(options_);
-[err] = ms_sbvar_command_line(perform_simulation);
-mexErrCheck('ms_sbvar_command_line simulate',err);
+disp('Simulation');
+options_ = set_ms_estimation_flags_for_other_mex(options_);
+oo_ = set_oo_w_estimation_output(options_, oo_);
 
-end
+% setup command line options
+opt = ['-simulate -seed ' num2str(options_.DynareRandomStreams.seed)];
+opt = [opt ' -ft ' options_.ms.output_file_tag];
+opt = [opt ' -ndraws ' num2str(options_.ms.mh_replic)];
+opt = [opt ' -burnin ' num2str(options_.ms.drop)];
+opt = [opt ' -thin ' num2str(options_.ms.thinning_factor)];
+opt = [opt ' -mh ' num2str(options_.ms.adaptive_mh_draws)];
 
-function opt=create_simulation_commandline(options_)
-
-opt = '-simulate';
-opt = [opt set_flag_if_exists(options_.ms, 'output_file_tag')];
-opt = [opt set_flag_if_exists(options_.ms, 'mh_replic')];
-opt = [opt set_flag_if_exists(options_.ms, 'drop')];
-opt = [opt set_flag_if_exists(options_.ms, 'thinning_factor')];
-opt = [opt set_flag_if_exists(options_.ms, 'adaptive_mh_draws')];
+% simulation
+[err] = ms_sbvar_command_line(opt);
+mexErrCheck('ms_sbvar_command_line simulate',err);
 
-if options_.DynareRandomStreams.seed
-    opt = [' -seed' num2str(options_.DynareRandomStreams.seed)];
-end
+options_.ms.simulation_file_tag = options_.ms.output_file_tag;
+options_ = initialize_ms_sbvar_options(M_, options_);
 end
diff --git a/matlab/ms-sbvar/ms_variance_decomposition.m b/matlab/ms-sbvar/ms_variance_decomposition.m
new file mode 100644
index 0000000000000000000000000000000000000000..556490ed5ddfeffd1ae21498c0d9de971a2ab37a
--- /dev/null
+++ b/matlab/ms-sbvar/ms_variance_decomposition.m
@@ -0,0 +1,67 @@
+function [options_, oo_]=ms_variance_decomposition(M_, options_, oo_)
+%function ms_variance_decomposition()
+% calls ms_variance_decomposition mex function
+%
+% INPUTS
+%    M_
+%    options_
+%    oo_
+%
+% OUTPUTS
+%    options_
+%    oo_
+%
+% SPECIAL REQUIREMENTS
+%    none
+
+% 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/>.
+
+disp('Variance Decomposition');
+options_ = set_ms_estimation_flags_for_other_mex(options_);
+options_ = set_ms_simulation_flags_for_other_mex(options_);
+oo_ = set_oo_w_estimation_output(options_, oo_);
+
+% NOTICE THAT VARIANCE DECOMPOSITION DEFAULTS TO USING THE MEAN, NOT MEDIAN OR BANDED
+
+opt = {options_.ms.output_file_tag, ...
+    'seed', options_.DynareRandomStreams.seed, ...
+    'horizon', options_.ms.horizon, ...
+    'filtered', options_.ms.filtered_probabilities, ...
+    'error_bands', options_.ms.error_bands, ...
+    'percentiles', options_.ms.error_band_percentiles, ...
+    'thin', options_.ms.thinning_factor, ...
+    'mean'};
+
+[err, vd] = mex_ms_variance_decomposition(opt{:},'free_parameters',oo_.ms.maxparams,'shocks',options_.ms.shock_draws);
+mexErrCheck('mex_ms_variance_decomposition ergodic ', err);
+plot_ms_variance_decomposition(M_,vd, 'Ergodic Variance Decomposition');
+
+[err, regime_vd] = mex_ms_variance_decomposition(opt{:},'free_parameters',oo_.ms.maxparams,'shocks',options_.ms.shock_draws,'regimes');
+mexErrCheck('mex_ms_variance_decomposition ergodic regimes', err);
+save([M_.fname '/' options_.ms.output_file_tag '_ergodic_vd.mat'], 'vd', 'regime_vd');
+
+if exist(options_.ms.load_mh_file,'file') > 0
+    [err, vd] = mex_ms_variance_decomposition(opt{:},'simulation_file',options_.ms.load_mh_file,'shocks',options_.ms.shocks_per_parameter,'parameter_uncertainty');
+    mexErrCheck('mex_ms_variance_decomposition bayesian ', err);
+
+    [err, regime_vd] = mex_ms_variance_decomposition(opt{:},'simulation_file',options_.ms.load_mh_file,'shocks',options_.ms.shocks_per_parameter,'parameter_uncertainty','regimes');
+    mexErrCheck('mex_ms_variance_decomposition bayesian regimes ', err);
+    save([M_.fname '/' options_.ms.output_file_tag '_bayesian_vd.mat'], 'vd', 'regime_vd');
+end
+options_ = initialize_ms_sbvar_options(M_, options_);
+end
diff --git a/matlab/ms-sbvar/plot_ms_forecast.m b/matlab/ms-sbvar/plot_ms_forecast.m
index 0ec63a7c3f0167c2f5928beea61419675d1fc041..2f7f7b352469cd02aeb4b84a1801336d74dfa253 100644
--- a/matlab/ms-sbvar/plot_ms_forecast.m
+++ b/matlab/ms-sbvar/plot_ms_forecast.m
@@ -1,8 +1,9 @@
-function plot_ms_forecast(forecast,title_)
+function plot_ms_forecast(M_,forecast,title_)
 % function [] = plot_ms_forecast(forecast,names)
 % plots the forecast from the output from a ms-sbvar
 %
 % INPUTS
+%   M_
 %   forecast should be in the form (percentile x horizon x nvar ), if banded otherwise
 %     ( horizon x nvar )
 %
@@ -26,9 +27,6 @@ function plot_ms_forecast(forecast,title_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
   
-  
-  global M_ oo_ options_
-
   nc = 2;
   nr = 2;
   nvars = M_.endo_nbr;
@@ -139,4 +137,4 @@ function plot_ms_forecast(forecast,title_)
   end
   
 
-end
\ No newline at end of file
+end
diff --git a/matlab/ms-sbvar/plot_ms_irf.m b/matlab/ms-sbvar/plot_ms_irf.m
index 2d37b2e99e24924d7fc4857eff63d575d3fd2140..0675b0966e3166cd17b1db775ec6598f7f0bcdb4 100644
--- a/matlab/ms-sbvar/plot_ms_irf.m
+++ b/matlab/ms-sbvar/plot_ms_irf.m
@@ -1,8 +1,9 @@
-function plot_ms_irf(irf,names,title_)
+function plot_ms_irf(M_,irf,names,title_)
 % function [] = plot_ms_irf(irf,names)
 % plots the impulse responses from the output from a ms-sbvar
 %
 % INPUTS
+%   M_
 %   irf should be in the form (percentile x horizon x (nvar x nvar)), if banded otherwise
 %     ( horizon x (nvar x nvar) )
 %
@@ -30,8 +31,6 @@ function plot_ms_irf(irf,names,title_)
 %
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
-  
-  global M_ oo_ options_
 
   nc = 2;
   nr = 2;
@@ -128,4 +127,4 @@ function plot_ms_irf(irf,names,title_)
     %suptitle(title_);
   end
   
-end
\ No newline at end of file
+end
diff --git a/matlab/ms-sbvar/plot_ms_probabilities.m b/matlab/ms-sbvar/plot_ms_probabilities.m
index 1c81b57d95edfe23812e2e38f88b4796678cd2f6..00694caaa9f053772ce782fe28e9e2367c0be0d8 100644
--- a/matlab/ms-sbvar/plot_ms_probabilities.m
+++ b/matlab/ms-sbvar/plot_ms_probabilities.m
@@ -1,4 +1,4 @@
-function plot_ms_probabilities(computed_probabilities)
+function plot_ms_probabilities(options_, computed_probabilities)
 % function plot_ms_probabilities(computed_probabilities)
 % Plots the regime probablities for each graph
 %
@@ -27,8 +27,6 @@ function plot_ms_probabilities(computed_probabilities)
 %
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
-
-  global M_ oo_ options_
   
   [T,num_grand_regimes] = size(computed_probabilities);
   num_chains = length(options_.ms.ms_chain);
@@ -65,4 +63,4 @@ function [chains] = iterate_chain(probs, t, chains, chain, num_chains)
       chains = iterate_chain(p, t, chains, chain+1, num_chains);
     end
   end
-end
\ No newline at end of file
+end
diff --git a/matlab/ms-sbvar/plot_ms_variance_decomposition.m b/matlab/ms-sbvar/plot_ms_variance_decomposition.m
index e7c8657943ab0d79ec4bc6f2e4b38d64fc27686e..d04b43c49ae3ecac2010a1f12b4fb8993fbf5b4f 100644
--- a/matlab/ms-sbvar/plot_ms_variance_decomposition.m
+++ b/matlab/ms-sbvar/plot_ms_variance_decomposition.m
@@ -1,7 +1,8 @@
-function plot_ms_variance_decomposition(vd, title_, varargin)
+function plot_ms_variance_decomposition(M_, vd, title_, varargin)
 % plot the variance decomposition of shocks
 %
 % Inputs
+%       M_
 %		shocks: matrix of the individual shocks Tx(KxK)with J=number of shocks
 %
 % Optional Inputs
@@ -31,8 +32,7 @@ function plot_ms_variance_decomposition(vd, title_, varargin)
 %
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
-	
-	global M_ oo_ options_
+
   nvars = M_.endo_nbr;
   endo_names = M_.endo_names;
 
@@ -145,4 +145,4 @@ function plot_ms_variance_decomposition(vd, title_, varargin)
 		%suptitle(title_);
 	end
 
-end
\ No newline at end of file
+end
diff --git a/matlab/ms-sbvar/set_flag_if_exists.m b/matlab/ms-sbvar/set_flag_if_exists.m
deleted file mode 100644
index cf6f2e24fba39791c9b9cb437b4b1c6e919e1bca..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/set_flag_if_exists.m
+++ /dev/null
@@ -1,97 +0,0 @@
-function flagstr = set_flag_if_exists(msoptions, flag)
-%function set_flag_if_exists()
-% adds a flag to the command line if it exists in the msoptions structure
-%
-% INPUTS
-%    msoptions:  options_.ms structure
-%    flag:       string representing the flag to check
-%
-% OUTPUTS
-%    none
-%
-% SPECIAL REQUIREMENTS
-%    none
-
-% 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 ~isfield(msoptions, flag)
-    flagstr = '';
-    return;
-end
-
-% Convert from Dynare options to ms-sbvar c-code options here
-% to allow for clear programming in the rest of the Dynare
-% matlab files
-switch flag
-    case 'output_file_tag'              % general option
-        cFlag = 'ft';
-    case 'convergence_starting_value'   % estimation options
-        cFlag = 'cb';
-    case 'convergence_ending_value'
-        cFlag = 'ce';
-    case 'convergence_increment_value'
-        cFlag = 'ci';
-    case 'max_iterations_starting_value'
-        cFlag = 'ib';
-    case 'max_iterations_increment_value'
-        cFlag = 'ii';
-    case 'max_block_iterations'
-        cFlag = 'mb';
-    case 'max_repeated_optimization_runs'
-        cFlag = 'repeat_max';
-    case 'function_convergence_criterion'
-        cFlag = 'repeat_tol_obj';
-    case 'parameter_convergence_criterion'
-        cFlag = 'repeat_tol_params';
-    case 'number_of_large_perturbations'
-        cFlag = 'random';
-    case 'number_of_small_perturbations'
-        cFlag = 'random_small';
-    case 'number_of_posterior_draws_after_perturbation'
-        cFlag = 'random_small_ndraws';
-    case 'max_number_of_stages'
-        cFlag = 'random_max';
-    case 'random_function_convergence_criterion'
-        cFlag = 'random_tol_obj';
-    case 'random_parameter_convergence_criterion'
-        cFlag = 'random_tol_params';
-    case 'mh_replic'                    % simulation options
-        cFlag = 'ndraws';
-    case 'drop'
-        cFlag = 'burnin';
-    case 'thinning_factor'
-        cFlag = 'thin';
-    case 'adaptive_mh_draws'
-        cFlag = 'mh';
-    case 'load_mh_file'                 % mdd options
-        cFlag = 'pf';
-    case 'mdd_proposal_draws'
-        cFlag = 'mdd_proposal_draws';
-    case 'mdd_use_mean_center'
-        cFlag = 'use_mean';
-    case 'filtered_probabilities'       % probabilities options
-        cFlag = 'filtered';
-    case 'real_time_smoothed_probabilities'
-        cFlag = 'real_time_smoothed';
-    otherwise
-        cFlag = flag;
-end
-
-flagstr = [' -' cFlag ' ' num2str(eval(['msoptions.' flag]))];
-
-end
\ No newline at end of file
diff --git a/matlab/ms-sbvar/set_input_options.m b/matlab/ms-sbvar/set_input_options.m
new file mode 100644
index 0000000000000000000000000000000000000000..f8b7e21ea80875f03cd0b55f8e326ad8a2cfb243
--- /dev/null
+++ b/matlab/ms-sbvar/set_input_options.m
@@ -0,0 +1,57 @@
+function options_=set_input_options(options_)
+%function set_input_options()
+% get mode_file from filename if it's not in mode_file
+%
+% INPUTS
+%    options_
+%
+% OUTPUTS
+%    options_
+%
+% SPECIAL REQUIREMENTS
+%    none
+
+% 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/>.
+
+% init file
+if ~isfield(options_.ms, 'init_file')
+    options_.ms.init_file = ['init_' options_.ms.output_file_tag '.dat'];
+end
+
+if ~exist(options_.ms.init_file, 'file')
+    error(['ERROR: cannot find init_file: ' options_.ms.init_file]);
+end
+
+% free parameter file
+if ~isfield(options_.ms, 'free_param_file')
+    options_.ms.free_param_file = ['est_free_' options_.ms.output_file_tag '.out'];
+end
+
+if ~exist(options_.ms.free_param_file, 'file')
+    error(['ERROR: cannot find free_param_file: ' options_.ms.free_param_file]);
+end
+
+% simulation file
+if ~isfield(options_.ms, 'load_mh_file')
+    options_.ms.load_mh_file = ['simulation_' options_.ms.output_file_tag '.out'];
+end
+
+if ~exist(options_.ms.load_mh_file, 'file')
+    error(['ERROR: Could not find Metropolis Hastings file: ' options_.ms.load_mh_file]);
+end
+end
diff --git a/matlab/ms-sbvar/set_ms_estimation_flags_for_other_mex.m b/matlab/ms-sbvar/set_ms_estimation_flags_for_other_mex.m
new file mode 100644
index 0000000000000000000000000000000000000000..0d179e2bc554be2f29dd047a2730d0de5f6057bd
--- /dev/null
+++ b/matlab/ms-sbvar/set_ms_estimation_flags_for_other_mex.m
@@ -0,0 +1,44 @@
+function options_=set_ms_estimation_flags_for_other_mex(options_)
+%function set_ms_estimation_flags_for_other_mex()
+% MS Sbvar Estimation
+%
+% INPUTS
+%    options_
+%
+% OUTPUTS
+%    options_
+%
+% SPECIAL REQUIREMENTS
+%    none
+
+% 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 ~isfield(options_.ms,'estimation_file_tag')
+    options_.ms.estimation_file_tag = options_.ms.output_file_tag;
+end
+options_.ms.free_param_file = ['est_free_' options_.ms.estimation_file_tag '.out'];
+options_.ms.init_file = ['init_' options_.ms.estimation_file_tag '.dat'];
+
+if ~exist(options_.ms.init_file,'file')
+    error(['ERROR: Could not find initialization file: ' options_.ms.init_file]);
+end
+
+if ~exist(options_.ms.free_param_file,'file')
+    error(['ERROR: Could not find free parameter file: ' options_.ms.free_param_file]);
+end
+end
diff --git a/matlab/ms-sbvar/set_ms_simulation_flags_for_other_mex.m b/matlab/ms-sbvar/set_ms_simulation_flags_for_other_mex.m
new file mode 100644
index 0000000000000000000000000000000000000000..a37da2cfa4ccfe5ebba7af352372a3b2fc037737
--- /dev/null
+++ b/matlab/ms-sbvar/set_ms_simulation_flags_for_other_mex.m
@@ -0,0 +1,39 @@
+function options_=set_ms_simulation_flags_for_other_mex(options_)
+%function set_ms_simulation_flags_for_other_mex()
+% MS Sbvar Estimation
+%
+% INPUTS
+%    options_
+%
+% OUTPUTS
+%    options_
+%
+% SPECIAL REQUIREMENTS
+%    none
+
+% 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 ~isfield(options_.ms,'simulation_file_tag')
+    options_.ms.simulation_file_tag = options_.ms.output_file_tag;
+end
+options_.ms.load_mh_file = ['simulation_' options_.ms.simulation_file_tag '.out'];
+
+if ~exist(options_.ms.load_mh_file,'file')
+    error(['ERROR: Could not find Metropolis Hastings file: ' options_.ms.load_mh_file]);
+end
+end
diff --git a/matlab/ms-sbvar/set_oo_w_estimation_output.m b/matlab/ms-sbvar/set_oo_w_estimation_output.m
new file mode 100644
index 0000000000000000000000000000000000000000..3bfd80125ebf5e5e6d0c3a25a248565d7ca5ffc1
--- /dev/null
+++ b/matlab/ms-sbvar/set_oo_w_estimation_output.m
@@ -0,0 +1,42 @@
+function oo_=set_oo_w_estimation_output(options_, oo_)
+%function set_oo_w_estimation_output()
+% places estimation output in oo_ structure
+%
+% INPUTS
+%    options_
+%    oo_
+%
+% OUTPUTS
+%    oo_
+%
+% SPECIAL REQUIREMENTS
+%    none
+
+% 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/>.
+
+
+oo_.ms.maxparams = load(options_.ms.free_param_file);
+oo_.ms.maxparams = oo_.ms.maxparams(3:end)';
+
+if ~isfield(oo_.ms, 'A0') || ~isfield(oo_.ms, 'Aplus') ...
+        || ~isfield(oo_.ms, 'Zeta') || ~isfield(oo_.ms, 'Q')
+    [err, oo_.ms.A0, oo_.ms.Aplus, oo_.ms.Zeta, oo_.ms.Q] = ...
+        mex_ms_convert_free_parameters(options_.ms.output_file_tag, oo_.ms.maxparams);
+    mexErrCheck('mex_ms_convert_free_parameters', err);
+end
+end
diff --git a/matlab/ms-sbvar/sz_prd.m b/matlab/ms-sbvar/sz_prd.m
deleted file mode 100644
index fb779068c6f6c2e55ea6f988db9ed0162fef6dcf..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/sz_prd.m
+++ /dev/null
@@ -1,800 +0,0 @@
-function sz_prd(M,options_)
-
-  global oo_;
-  %==========================================================================
-  %== Directory structure
-  %==========================================================================
-
-  %generation of mhm file
-  %generateMHM_input(options_);
-
-  ms_root = strrep(which('ms_sbvar'),'/ms_sbvar.m','');
-
-  % path for C executables
-  %c_path='./c-executables';
-  c_path = [ms_root '../../ms-sbvar'];
-
-  % path for Markov specification
-  %m_spec_path=[ms_root '/switching_specification'];
-
-  % path for MHM specification
-  %mhm_spec_path=[ms_root '/mhm_specification'];
-
-  %==========================================================================
-  %== Processing control
-  %==========================================================================
-  % 1 to use standalone, 0 to use mex
-  options_.ms.standalone = 0;
-
-  %==========================================================================
-  %== Output control
-  %==========================================================================
-  % tag for output files %create an option
-  %options_.ms.output_file_tag='test_2v';
-
-  % 1 to create init_<options_.ms.output_file_tag>.dat, 0 otherwise
-  %options_.ms.create_initialization_file = 1; %0 originally
-
-  % 1 to perform estimation, 0 otherwise
-  %options_.ms.estimate_msmodel = 1;
-
-  % 1 to perform estimation, 0 otherwise
-  %options_.ms.compute_mdd = 1;
-
-  % 1 to compute probabilites, 0 otherwise
-  %options_.ms.compute_probabilities = 1;%1 in the original
-
-  % 1 to Prints draws of the posterior
-  %options_.ms.print_draws = 1;
-  %options_.ms.n_draws=1000;
-  %options_.ms.thinning_factor=1;
-
-  %==========================================================================
-  %== Markov Process Specification File
-  %==========================================================================
-  %options_.ms.markov_file = 'specification_2v2c.dat'; %create an option
-  markov_file = [options_.ms.markov_file '.dat'];
-
-  %==========================================================================
-  %== Markov Process Specification File
-  %==========================================================================
-  %options_.ms.mhm_file = 'MHM_input.dat';
-
-
-  %mhm_file = '/MHM_input.dat';
-  %options_.ms.proposal_draws = 100000;
-
-  %==========================================================================
-  %== Var Specification
-  %==========================================================================
-  % Number of options_.ms.nlags
-  %options_.ms.nlags = 4;
-
-  % Var restriction function
-  %options_.ms.restriction_fname = 'ftd_upperchol3v'; %create an option
-
-
-  %==========================================================================
-  %== BVAR prior
-  %==========================================================================
-
-  %=== The following mu is effective only if indxPrior==1.
-  mu = zeros(6,1);   % hyperparameters
-  if length(options_.ms.coefficients_prior_hyperparameters) ~= 6
-    error('When specifying the coefficients_prior_hyperparameters, you must pass a vector of 6 numbers')
-  end
-  mu = options_.ms.coefficients_prior_hyperparameters;
-  mu = reshape(mu,1,6);
-
-  % mu(1): overall tightness for A0 and Aplus
-  % mu(2): relative tightness for Aplus
-  % mu(3): relative tightness for the constant term
-  % mu(4): tightness on lag decay.  (1.2 - 1.5 faster decay produces better
-  % inflation forecasts
-  % mu(5): weight on nvar sums of coeffs dummy observations (unit roots).
-  % mu(6): weight on single dummy initial observation including constant
-  % (cointegration, unit roots, and stationarity).
-
-  % Alpha on p. 66 for squared time-varying structural shock lambda.
-  galp = 1.0; 
-
-  % Beta on p. 66 for squared time-varying structural shock lambda.
-  gbeta = 1.0;   
-
-  % Case 3 (no state change across options_.ms.nlags (l) but allows all variables for a
-  % given lag to switch states). Normal prior variance for glamda
-  % (nvar-by-nvar for each state) for different variables in lagged D+.  See
-  % p.71v.  
-  gsig2_lmdm = 50^2;  
-
-
-  %==========================================================================
-  %== Data
-  %==========================================================================
-  % Read in data to produce rectangular array named xdd.  Each column is one
-  % data series.
-  %load ./data/artificial_data
-  xdd=options_.data;
-
-  % Information about timing of the data for consistancy checks
-  % quarters (4) or months (12)
-  %q_m = 4;  
-  %options_.ms.freq = 4;
-  q_m = options_.ms.freq;
-  % beginning year in data set
-  %yrBin=1978;   
-  %options_.ms.initial_year = 1959;
-  yrBin=options_.ms.initial_year;
-  % beginning quarter or month in data set
-  %qmBin=1;  
-  %options_.ms.initial_subperiod = 1;
-  qmBin=options_.ms.initial_subperiod;
-  % final year in data set
-  %yrFin=2007;  
-  %options_.ms.final_year = 2005;
-  yrFin=options_.ms.final_year;
-  % final month or quarter in data set
-  %qmFin=4;  
-  %options_.ms.final_subperiod = 4;
-  qmFin=options_.ms.final_subperiod;
-  % first year to use in estimation
-  %yrStart=yrBin;
-  yrStart=options_.ms.initial_year;
-  % first quarter or month to use in estimation
-  %qmStart=qmBin;
-  qmStart=options_.ms.initial_subperiod;
-  % last year to use in estimation
-  %yrEnd=yrFin;
-  yrEnd=options_.ms.final_year;
-  % last quater or month to use in estimation
-  %qmEnd=qmFin;
-  qmEnd=options_.ms.final_subperiod;
-  % Log variables in xdd
-  logindx = []; 
-
-  % Convert percent to decimal in xdd
-  pctindx = [];  
-
-  % Select the variable to use and rearrange columns if desired
-  %vlist = [3 1 2];
-  %options_.ms.vlist = [1 2 3];
-  options_.ms.vlist = [1:size(options_.varobs,1)];
-  vlist1=options_.ms.vlist;
-
-  %==========================================================================
-  %== Linux or Windows system - differences in some naming conventions
-  %==========================================================================
-  use_linux = 1;
-
-
-  %==========================================================================
-  %== Random number seed for selecting starting point in constant parameter 
-  %== optimization.
-  %==========================================================================
-  % Set to zero to set from clock
-  %seednumber = 7910;
-  %if seednumber
-  %    randn('state',seednumber);
-  %    rand('state',seednumber);
-  %else
-  %    rand('state',fix(100*sum(clock)));
-  %    randn('state',1000000000*rand);
-  %    rand('state',1000000000*rand);
-  %end
-  seednumber = options_.DynareRandomStreams.seed;
-
-  %==========================================================================
-  %==========================================================================
-  %==========================================================================
-  %== Beginning of code.  Modify below at own risk.
-  %==========================================================================
-
-  % options that may at some point become user specified
-  %indxC0Pres = 0;    % 1: cross-A0-and-A+ restrictions; 0: idfile_const is all we have
-  indxC0Pres =options_.ms.cross_restrictions;
-  % Example for indxOres==1: restrictions of the form P(t) = P(t-1).
-  %Rform = 0;         % 1: contemporaneous recursive reduced form; 0: restricted (non-recursive) form
-  Rform =options_.ms.contemp_reduced_form;
-  % % % Pseudo = 0;  % 1: Pseudo forecasts; 0: real time forecasts
-  %indxPrior = 1;     % 1: Bayesian prior; 0: no prior
-  indxPrior =options_.ms.bayesian_prior;
-  %indxDummy = indxPrior;  % 1: add dummy observations to the data; 0: no dummy added.
-  indxDummy = options_.ms.bayesian_prior;
-  %ndobs = 0;         % No dummy observations for xtx, phi, fss, xdatae, etc.  Dummy observations are used as an explicit prior in fn_rnrprior_covres_dobs.m.
-  %ndobs =options_.ms.dummy_obs;
-  %if indxDummy
-  %   ndobs=nvar+1;         % number of dummy observations
-  %else
-  %   ndobs=0;    % no dummy observations
-  %end
-
-  % 
-  hpmsmd = [0.0; 0.0];
-  indxmsmdeqn = [0; 0; 0; 0];  %This option disenable using this in fn_rnrprior_covres_dobs.m
-
-  nStates = -1;
-
-
-
-
-  %==========================================================================
-  %== Create initialization file
-  %==========================================================================
-  if options_.ms.create_initialization_file == 1
-    %======================================================================
-    %== Check and setup data
-    %======================================================================
-    % log data
-    xdd(:,logindx) = log(xdd(:,logindx));
-
-    % convert percentage to decimal
-    xdd(:,pctindx)=.01*xdd(:,pctindx);
-
-    if (q_m ~= 12) && (q_m ~= 4)
-      disp('Warning: data must be monthly or quarterly!')
-      return
-    end
-
-    % number of data points
-    nData=(yrFin-yrBin)*q_m + (qmFin-qmBin+1);
-    % number of data points in estimation sample
-    nSample=(yrEnd-yrStart)*q_m + (qmEnd-qmEnd+1);
-    % number of periods not used at beginning of data (non-negative number)
-    nStart=(yrStart-yrBin)*q_m + (qmStart-qmBin);
-    % number of periods not used at end of data (non-positive number)
-    nEnd=(yrEnd-yrFin)*q_m + (qmEnd-qmFin);
-
-    if (nEnd > 0) || (nStart < 0)
-      disp('Warning: desired estimation period not in data set!')
-      return
-    end
-    if (nSample <= 0)
-      disp('Warning: no data points in estimation period!')
-      return
-    end
-
-    % reorder variables and create estimation data set
-    xdgel=xdd(nStart+1:nData+nEnd,vlist1);
-
-    % bad data points
-    baddata = find(isnan(xdgel));
-    if ~isempty(baddata)
-      disp('Warning: some data for estimation period are unavailable.')
-      return
-    end
-
-    % set nvar and nexo
-    nvar=size(xdgel,2);
-    nexo=1;
-
-    % Arranged data information, WITHOUT dummy obs when 0 after mu is used.
-    % See fn_rnrprior_covres_dobs.m for using the dummy observations as part of
-    % an explicit prior.
-    [xtx,xty,yty,fss,phi,y,ncoef,xr,Bh] = fn_dataxy(nvar,options_.ms.nlags,xdgel,mu,0,nexo);
-
-
-    %======================================================================
-    %== Linear Restrictions
-    %======================================================================
-    if Rform
-      Ui=cell(nvar,1); Vi=cell(ncoef,1);
-      for kj=1:nvar
-        Ui{kj} = eye(nvar);  Vi{kj} = eye(ncoef);
-      end
-    else
-      [Ui,Vi,n0,np,ixmC0Pres] = feval(options_.ms.restriction_fname,nvar,nexo,options_.ms);
-      if min(n0)==0
-        disp('A0: restrictions give no free parameters in one of equations')
-        return
-      elseif min(np)==0
-        disp('Aplus: Restrictions in give no free parameters in one of equations')
-        return
-      end
-    end
-
-
-    %======================================================================
-    %== Estimation
-    %======================================================================
-    if indxPrior
-      %*** Obtains asymmetric prior (with no linear restrictions) with dummy observations as part of an explicit prior (i.e,
-      %      reflected in Hpmulti and Hpinvmulti).  See Forecast II, pp.69a-69b for details.
-      if 1  % Liquidity effect prior on both MS and MD equations.
-        [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] = fn_rnrprior_covres_dobs(nvar,q_m,options_.ms.nlags,xdgel,mu,indxDummy,hpmsmd,indxmsmdeqn);
-      else
-        [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] = fn_rnrprior(nvar,q_m,options_.ms.nlags,xdgel,mu);
-      end
-
-      %*** Combines asymmetric prior with linear restrictions
-      [Ptld,H0invtld,Hpinvtld] = fn_rlrprior(Ui,Vi,Pi,H0multi,Hpmulti,nvar);
-
-      %*** Obtains the posterior matrices for estimation and inference
-      [Pmat,H0inv,Hpinv] = fn_rlrpostr(xtx,xty,yty,Ptld,H0invtld,Hpinvtld,Ui,Vi);
-    else
-      %*** Obtain the posterior matrices for estimation and inference
-      [Pmat,H0inv,Hpinv] = fn_dlrpostr(xtx,xty,yty,Ui,Vi);
-    end
-
-    if Rform
-      %*** Obtain the ML estimate
-      A0hatinv = chol(H0inv{1}/fss);   % upper triangular but lower triangular choleski
-      A0hat=inv(A0hatinv);
-
-      Aphat = Pmat{1}*A0hat;
-    else
-      %*** Obtain the ML estimate
-      %   load idenml
-      x = 10*rand(sum(n0),1);
-      H0 = eye(sum(n0));
-      crit = 1.0e-9;
-      nit = 10000;
-      %
-      tic
-      [fhat,xhat,grad,Hhat,itct,fcount,retcodehat] = csminwel('fn_a0freefun',x,H0,'fn_a0freegrad',crit,nit,Ui,nvar,n0,fss,H0inv);
-    endtime = toc;
-
-    A0hat = fn_tran_b2a(xhat,Ui,nvar,n0);
-
-    xhat = fn_tran_a2b(A0hat,Ui,nvar,n0);
-    [Aphat,ghat] = fn_gfmean(xhat,Pmat,Vi,nvar,ncoef,n0,np);
-    if indxC0Pres
-      Fhatur0P = Fhat;  % ur: unrestriced across A0 and A+
-      for ki = 1:size(ixmC0Pres,1)   % loop through the number of equations in which
-        % cross-A0-A+ restrictions occur. See St. Louis Note p.5.
-        ixeq = ixmC0Pres{ki}(1,1);   % index for the jth equation in consideration.
-        Lit = Vi{ixeq}(ixmC0Pres{ki}(:,2),:);  % transposed restriction matrix Li
-        % V_j(i,:) in f_j(i) = V_j(i,:)*g_j
-        ci = ixmC0Pres{ki}(:,4) .* A0hat(ixmC0Pres{ki}(:,3),ixeq);
-        % s * a_j(h) in the restriction f_j(i) = s * a_j(h).
-        LtH = Lit/Hpinv{ixeq};
-        HLV = LtH'/(LtH*Lit');
-        gihat = Vi{ixeq}'*Fhatur0P(:,ixeq);
-        Aphat(:,ixeq) = Vi{ixeq}*(gihat + HLV*(ci-Lit*gihat));
-      end
-    end
-  end
-
-
-  %======================================================================
-  %== Create matlab initialization file
-  %======================================================================
-  matlab_filename = ['matlab_',options_.ms.output_file_tag,'.prn'];
-  fidForC = fopen(matlab_filename,'w');
-
-  fprintf(fidForC,'\n%s\n','//== gxia: alpha parameter for gamma prior of xi ==//');
-  fprintf(fidForC,' %20.15f ', galp);
-  fprintf(fidForC, '\n\n');
-
-  fprintf(fidForC,'\n%s\n','//== gxib: beta parameter for gamma prior of xi ==//');
-  fprintf(fidForC,' %20.15f ', gbeta);
-  fprintf(fidForC, '\n\n');
-
-  fprintf(fidForC,'\n%s\n','//== glamdasig: sigma parameter for normal prior of lamda ==//');
-  fprintf(fidForC,' %20.15f ', sqrt(gsig2_lmdm));
-  fprintf(fidForC, '\n\n');
-
-  %=== lags, nvar, nStates, sample size (excluding options_.ms.nlags where, with dummyies, fss=nSample-options_.ms.nlags+ndobs).
-  fprintf(fidForC,'\n%s\n','//== lags, nvar, nStates, T ==//');
-  fprintf(fidForC,' %d  %d  %d  %d\n\n\n',options_.ms.nlags, nvar, nStates, fss);
-
-  %=== A0hat nvar-by-nvar from the constant VAR.
-  fprintf(fidForC,'\n%s\n','//== A0hat: nvar-by-nvar ==//');
-  indxFloat = 1;
-  xM = A0hat;
-  nrows = nvar;
-  ncols = nvar;
-  fn_fprintmatrix(fidForC, xM, nrows, ncols, indxFloat)
-
-  %=== Aphat ncoef-by-nvar from the constant VAR.
-  %=== Each column of Aphat is in the order of [nvar variables for 1st lag, ..., nvar variables for last lag, constant term].
-  fprintf(fidForC,'\n%s\n','//== Aphat: ncoef(lags*nvar+1)-by-nvar ==//');
-  indxFloat = 1;
-  xM = Aphat;
-  nrows = ncoef;
-  ncols = nvar;
-  fn_fprintmatrix(fidForC, xM, nrows, ncols, indxFloat)
-
-  %=== n0const: nvar-by-1, whose ith element represents the number of free A0 parameters in ith equation for the case of constant parameters.
-  fprintf(fidForC,'\n%s\n','//== n0const: nvar-by-1 ==//');
-  indxFloat = 0;
-  xM = n0;
-  nrows = 1;
-  ncols = nvar;
-  fn_fprintmatrix(fidForC, xM', nrows, ncols, indxFloat)
-
-  %=== npconst: nvar-by-1, whose ith element represents the number of free A+ parameters in ith equation for the case of constant parameters.
-  fprintf(fidForC,'\n%s\n','//== npconst: nvar-by-1 ==//');
-  indxFloat = 0;
-  xM = np;
-  nrows = 1;
-  ncols = nvar;
-  fn_fprintmatrix(fidForC, xM', nrows, ncols, indxFloat)
-
-  %=== Uiconst: nvar-by-1 cell.  In each cell, nvar-by-qi orthonormal basis for the null of the ith
-  %           equation contemporaneous restriction matrix where qi is the number of free parameters.
-  %           With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector
-  %           of total original parameters and bi is a vector of free parameters. When no
-  %           restrictions are imposed, we have Ui = I.  There must be at least one free
-  %           parameter left for the ith equation.
-  fprintf(fidForC,'\n%s\n','//== Uiconst: cell(nvar,1) and nvar-by-n0const(i) for the ith cell (equation) ==//');
-  for i_=1:nvar
-    fn_fprintmatrix(fidForC, Ui{i_}, nvar, n0(i_), 1);
-  end
-
-  %=== Viconst: nvar-by-1 cell.  In each cell, k-by-ri orthonormal basis for the null of the ith
-  %           equation lagged restriction matrix where k is a total of exogenous variables and
-  %           ri is the number of free parameters. With this transformation, we have fi = Vi*gi
-  %           or Vi'*fi = gi where fi is a vector of total original parameters and gi is a
-  %           vector of free parameters. There must be at least one free parameter left for
-  %           the ith equation.
-  fprintf(fidForC,'\n%s\n','//== Viconst: cell(nvar,1) and ncoef-by-n0const(i) for the ith cell (equation) ==//');
-  for i_=1:nvar
-    fn_fprintmatrix(fidForC, Vi{i_}, ncoef, np(i_), 1);
-  end
-
-  %=== H0barconstcell: cell(nvar,1) (equations) and n-by-n for each cell (equaiton).
-  %=== H0barconst:  prior covariance matrix for each column of A0 under asymmetric prior (including SZ dummy obs.) with NO linear restrictions imposed yet.
-  fprintf(fidForC,'\n%s\n','//== H0barconstcell: cell(nvar,1) and n-by-n for the ith cell (equation) ==//');
-  for i_=1:nvar
-    fn_fprintmatrix(fidForC, H0multi(:,:,i_), nvar, nvar, 1);
-  end
-
-  %=== Hpbarconstcell: cell(nvar,1) (equations) and ncoef-by-ncoef for each cell (equaiton).
-  %=== Hpbarconst:  prior covariance matrix for each column of A+ under asymmetric prior (including SZ dummy obs.) with NO linear restrictions imposed yet.
-  fprintf(fidForC,'\n%s\n','//== Hpbarconstcell: cell(nvar,1) and ncoef-by-ncoef for the ith cell (equation) ==//');
-  for i_=1:nvar
-    fn_fprintmatrix(fidForC, Hpmulti(:,:,i_), ncoef, ncoef, 1);
-  end
-
-  %=== phi:  X; T-by-k; column: [nvar for 1st lag, ..., nvar for last lag, other exogenous terms, const term]
-  fprintf(fidForC,'\n%s\n','//== Xright -- X: T-by-ncoef ==//');
-  xM = phi;
-  nrows = fss;
-  ncols = ncoef;
-  for ki=1:nrows
-    for kj=1:ncols
-      fprintf(fidForC,' %20.15f ',xM((kj-1)*nrows+ki));
-      if (kj==ncols)
-        fprintf(fidForC,'\n');
-      end
-    end
-    if (ki==nrows)
-      fprintf(fidForC,'\n\n');
-    end
-  end
-
-  %=== y:    Y: T-by-nvar where T=fss
-  fprintf(fidForC,'\n%s\n','//== Yleft -- Y: T-by-nvar ==//');
-  xM = y;
-  nrows = fss;
-  ncols = nvar;
-  for ki=1:nrows
-    for kj=1:ncols
-      fprintf(fidForC,' %20.15f ',xM((kj-1)*nrows+ki));
-      if (kj==ncols)
-        fprintf(fidForC,'\n');
-      end
-    end
-    if (ki==nrows)
-      fprintf(fidForC,'\n\n');
-    end
-  end
-
-  fclose(fidForC);
-
-  %======================================================================
-  %== Create C initialization filename
-  %======================================================================
-  ms_write_markov_file(markov_file,M,options_)
-  if options_.ms.standalone == 1
-    if use_linux == 1
-      create_init_file=[c_path,'/sbvar_init_file ',matlab_filename,' ',markov_file,' ',options_.ms.output_file_tag];
-      system(create_init_file); %Run operating system command and return result
-    else
-      create_init_file=[c_path,'\sbvar_init_file.exe ',matlab_filename,' ',markov_file,' ',options_.ms.output_file_tag];
-      dos(create_init_file)            
-    end
-  else           
-    create_init_file=[matlab_filename,' ',markov_file,' ',options_.ms.output_file_tag];
-    [err] = ms_sbvar_create_init_file(create_init_file);
-    mexErrCheck('ms_sbvar_create_init_file',err);
-  end
-end
-
-
-
-%==========================================================================
-%== Perform estimation
-%==========================================================================
-if options_.ms.mode_compute > 0
-  opt = ['-estimate -seed ' num2str(seednumber)];
-  opt = [opt ' -random ' num2str(options_.ms.estimate.random)];
-  opt = [opt ' -random_max ' num2str(options_.ms.estimate.random_max)];
-  opt = [opt ' -random_tol_obj ' num2str(options_.ms.estimate.random_tol_obj)];
-  opt = [opt ' -random_tol_parms ' num2str(options_.ms.estimate.random_tol_parms)];
-  opt = [opt ' -cb ' num2str(options_.ms.estimate.cb)];
-  opt = [opt ' -ce ' num2str(options_.ms.estimate.ce)];
-  opt = [opt ' -ci ' num2str(options_.ms.estimate.ci)];
-  opt = [opt ' -ib ' num2str(options_.ms.estimate.ib)];
-  opt = [opt ' -ii ' num2str(options_.ms.estimate.ii)];
-  opt = [opt ' -mb ' num2str(options_.ms.estimate.mb)];
-  if ischar(options_.ms.mode_file) > 0 && exist(options_.ms.mode_file) > 0
-    opt = [opt ' -fp ' options_.ms.mode_file];
-  end
-  opt = [opt ' -ft '];
-
-  if options_.ms.standalone == 1
-    if use_linux == 1
-      perform_estimation=[c_path,['/sbvar_commandline ' opt ],options_.ms.output_file_tag];
-      system(perform_estimation);
-    else
-      perform_estimation=[c_path,['\sbvar_commandline.exe' opt],options_.ms.output_file_tag];
-      dos(perform_estimation)            
-    end
-  else
-    perform_estimation=[opt ,options_.ms.output_file_tag];
-    [err] = ms_sbvar_command_line(perform_estimation);
-    mexErrCheck('ms_sbvar_command_line estimation',err);
-  end
-end
-
-%==========================================================================
-%== Simulation
-%==========================================================================
-if options_.ms.mh_replic > 0
-  opt = ['-simulate -seed ' num2str(seednumber)];
-  opt = [opt ' -ndraws ' num2str(options_.ms.mh_replic)]; 
-  opt = [opt ' -burnin ' num2str(options_.ms.drop)];  
-  opt = [opt ' -mh ' num2str(options_.ms.adapative_mh_draws)];
-  opt = [opt ' -thin ' num2str(options_.ms.thinning_factor)];
-  if options_.ms.mode_compute == 0 && ischar(options_.ms.mode_file) > 0 && exist(options_.ms.mode_file) > 0
-    opt = [opt ' -fp ' options_.ms.mode_file];
-  end
-  opt = [opt ' -ft '];
-
-  if options_.ms.standalone == 1
-    if use_linux == 1
-      print_draws=[c_path,'/sbvar_commandline ', opt,options_.ms.output_file_tag];
-      system(print_draws);
-    else
-      print_draws=[c_path,'\sbvar_commandline.exe ',opt ,options_.ms.output_file_tag];
-      system(print_draws);
-    end
-  else
-    print_draws=[opt,options_.ms.output_file_tag];
-    [err] = ms_sbvar_command_line(print_draws);
-    mexErrCheck('ms_sbvar_command_line simulation',err);
-  end
-end
-
-%==========================================================================
-%== Compute marginal data density
-%==========================================================================
-if options_.ms.compute_mdd == 1
-  opt = ['-mdd -seed ' num2str(seednumber)];
-  opt = [opt ' -d ' num2str(options_.ms.mdd_proposal_draws)];
-  opt = [opt ' -pt ' num2str(options_.ms.mdd_proposal_type(1))];
-  opt = [opt ' -l ' num2str(options_.ms.mdd_proposal_type(2))];
-  opt = [opt ' -h ' num2str(options_.ms.mdd_proposal_type(3))];
-  % pf file?
-  if options_.ms.mdd_use_mean_center
-    opt = [opt ' -use_mean'];
-  end
-  opt = [opt ' -ft '];
-  if options_.ms.standalone == 1
-    if use_linux == 1
-      compute_mdd=[c_path,'/sbvar_commandline ',opt,options_.ms.output_file_tag];
-      system(compute_mdd);
-    else
-      compute_mdd=[c_path,'\sbvar_commandline.exe ',opt,options_.ms.output_file_tag];
-      system(compute_mdd);
-    end
-  else
-    compute_mdd=[opt,options_.ms.output_file_tag];
-    [err] = ms_sbvar_command_line(compute_mdd);
-    mexErrCheck('ms_sbvar_command_line marginal data density',err);
-  end
-
-  % grab the muller/bridge log mdd from the output file
-  mull_exp = 'Muller \w+\(\w+\) \= (\d+.\d+e\+\d+)';
-  bridge_exp = 'Bridge \w+\(\w+\) \= (\d+.\d+e\+\d+)';
-  bridge_mdd = -1; muller_mdd = -1;
-  mdd_filename = ['mdd_t' num2str(options_.ms.mdd_proposal_type(1)) '_' options_.ms.output_file_tag '.out'];
-  if exist(mdd_filename,'file')
-    mdd_fid = fopen(mdd_filename);
-    tline = fgetl(mdd_fid);
-    while ischar(tline)
-      mull_tok = regexp(tline,mull_exp,'tokens');
-      bridge_tok = regexp(tline,bridge_exp,'tokens');
-      if (~isempty(mull_tok))
-        muller_mdd = str2double(mull_tok{1}{1});
-      end
-      if (~isempty(bridge_tok))
-        bridge_mdd = str2double(bridge_tok{1}{1});
-      end
-      tline = fgetl(mdd_fid);
-    end
-    oo_.ms.mueller_log_mdd = muller_mdd;
-    oo_.ms.bridged_log_mdd = bridge_mdd;
-  end
-end
-
-
-
-%==========================================================================
-%== Compute posterior mode regime probabilities
-%==========================================================================
-if options_.ms.compute_probabilities == 1
-  opt = ['-probabilities -seed ' num2str(seednumber)];
-  if options_.ms.filtered_probabilities
-    opt = [opt ' -filtered' ];
-    prob_out_file = ['filtered_' options_.ms.output_file_tag '.out'];
-  else
-    opt = [opt ' -smoothed' ];
-    prob_out_file = ['smoothed_' options_.ms.output_file_tag '.out'];
-  end
-  if options_.ms.real_time_smoothed_probabilities
-    opt = [opt ' -real_time_smoothed' ];
-  end
-  opt = [opt ' -ft '];
-  if options_.ms.standalone == 1
-    if use_linux == 1
-      compute_prob=[c_path,'/sbvar_commandline ',opt,options_.ms.output_file_tag];
-      system(compute_prob);
-    else
-      compute_prob=[c_path,'\sbvar_commandline.exe ',opt,options_.ms.output_file_tag];
-      system(compute_prob);
-    end
-  else
-    compute_prob=[opt,options_.ms.output_file_tag];
-    [err] = ms_sbvar_command_line(compute_prob);
-    mexErrCheck('ms_sbvar_command_line probabilities',err);
-
-    % now we want to plot the probabilities for each chain
-    computed_probabilities = load(prob_out_file);
-    plot_ms_probabilities(computed_probabilities);
-  end
-end
-
-%==========================================================================
-%== Set the mode file and simulation file
-%==========================================================================
-
-if options_.ms.mode_compute == 0 && ischar(options_.ms.mode_file) > 0 && exist(options_.ms.mode_file,'file') > 0
-  mode_file = options_.ms.mode_file;
-else
-  mode_file = ['est_free_' options_.ms.output_file_tag '.out'];
-end
-if exist(mode_file,'file')
-  maxparams = load(mode_file);
-  maxparams = maxparams(3:end)';
-  [A0hat,Aplushat,Zetahat,Qhat] = mex_ms_convert_free_parameters(options_.ms.output_file_tag,maxparams);
-  oo_.ms.A0 = A0hat; oo_.ms.Aplus = Aplushat; oo_.ms.Zeta = Zetahat; oo_.ms.Q = Qhat;
-  mexErrCheck('mex_ms_convert_free_parameters',0);
-else
-  error('There is no mode file, can not continue');
-end
-
-if ischar(options_.ms.load_mh_file) > 0
-  if ~exist(options_.ms.load_mh_file,'file')
-    error('The Metropolis Hastings File supplied does not exist');
-  end
-  mh_file = options_.ms.load_mh_file;
-else
-  mh_file = ['simulation_' options_.ms.output_file_tag '.out'];  
-end
-
-fname = M.fname;
-if ~exist(fname, 'dir')
-  mkdir('.',fname);
-end
-
-%==========================================================================
-%== Compute Impulse Responses
-%==========================================================================
-
-if options_.ms.irf > 0
-  disp(['IRFs']); 
-
-  opt = {options_.ms.output_file_tag, 'seed', seednumber,'horizon',options_.ms.irf, 'percentiles', options_.ms.percentiles, ...
-  'filtered',options_.ms.irf_filtered, 'thin', options_.ms.irf_thinning_factor };
-
-  [irf] = mex_ms_irf(opt{:}, 'free_parameters', maxparams, 'shocks_per_parameter', options_.ms.irf_shock_draws);
-  if isscalar(irf)
-    mexErrCheck('mex_ms_irf ergodic ',irf);
-  end
-  plot_ms_irf(irf,options_.varobs,'Ergodic Impulse Responses');
-  [regime_irfs] = mex_ms_irf(opt{:}, 'free_parameters',maxparams,'shocks_per_parameter', options_.ms.irf_shock_draws,'regimes');
-  if isscalar(regime_irfs)
-    mexErrCheck('mex_ms_irf ergodic regimes ',regime_irfs);
-  end
-  save([fname '/' fname '_ergodic_irf.mat'], 'irf', 'regime_irfs');
-
-  if exist(mh_file,'file') > 0
-    [irf] = mex_ms_irf(opt{:}, 'shocks_per_parameter', options_.ms.irf_shocks_per_parameter, ...
-    'parameter_uncertainty','simulation_file',mh_file);
-    if isscalar(irf)
-      mexErrCheck('mex_ms_irf bayesian ',irf);
-    end
-    plot_ms_irf(irf,options_.varobs,'Impulse Responses w/ Parameter Uncertainty');
-    [regime_irfs] = mex_ms_irf(opt{:}, 'shocks_per_parameter', options_.ms.irf_shocks_per_parameter, ...
-    'simulation_file',mh_file,'parameter_uncertainty','regimes');
-    if isscalar(regime_irfs)
-      mexErrCheck('mex_ms_irf bayesian regimes ',regime_irfs);
-    end
-    save([fname '/' fname '_bayesian_irf.mat'], 'irf', 'regime_irfs');
-  end
-end
-
-%==========================================================================
-%== Compute Forecasts
-%==========================================================================
-if options_.ms.forecast > 0
-  disp('Forecasts'); 
-
-  opt = {options_.ms.output_file_tag, 'seed', seednumber,'horizon',options_.ms.forecast, 'percentiles', options_.ms.percentiles, ...
-  'thin', options_.ms.forecast_thinning_factor, 'data',options_.ms.forecast_data_obs };
-
-  [forecast] = mex_ms_forecast(opt{:},'free_parameters',maxparams,'shocks_per_parameter', options_.ms.forecast_shock_draws);
-  if isscalar(forecast)
-    mexErrCheck('mex_ms_forecast ergodic ',forecast);
-  end
-  plot_ms_forecast(forecast,'Forecast');
-  [regime_forecast] = mex_ms_forecast(opt{:},'free_parameters',maxparams,'shocks_per_parameter', options_.ms.forecast_shock_draws,'regimes');
-  if isscalar(regime_forecast)
-    mexErrCheck('mex_ms_forecast ergodic regimes',regime_forecast);
-  end
-  save([fname '/' fname '_ergodic_forecast.mat'], 'forecast', 'regime_forecast');  
-
-  if exist(mh_file,'file') > 0
-    [forecast] = mex_ms_forecast(opt{:},'free_parameters',maxparams,'shocks_per_parameter', options_.ms.forecast_shocks_per_parameter, ...
-    'simulation_file',mh_file,'parameter_uncertainty');
-    plot_ms_forecast(forecast,'Forecast w/ Parameter Uncertainty');
-    if isscalar(forecast)
-      mexErrCheck('mex_ms_forecast bayesian ',forecast);
-    end
-    [regime_forecast] = mex_ms_forecast(opt{:},'free_parameters',maxparams,'shocks_per_parameter', options_.ms.forecast_shocks_per_parameter, ...
-    'simulation_file',mh_file,'parameter_uncertainty','regimes');
-    if isscalar(regime_forecast)
-      mexErrCheck('mex_ms_forecast bayesian regimes ',regime_forecast);
-    end
-    save([fname '/' fname '_bayesian_forecast.mat'], 'forecast', 'regime_forecast'); 
-  end
-end
-
-%==========================================================================
-%== Compute Variance Decomposition
-%==========================================================================
-if options_.ms.variance_decomposition > 0
-  disp('Variance Decomposition');
-
-  % NOTICE THAT VARIANCE DECOMPOSITION DEFAULTS TO USING THE MEAN, NOT MEDIAN OR BANDED
-
-  opt = {options_.ms.output_file_tag, 'seed', seednumber,'horizon',options_.ms.forecast, ...
-  'thin', options_.ms.vd_thinning_factor,'filtered',options_.ms.vd_filtered, 'mean'};
-
-  [vd] = mex_ms_variance_decomposition(opt{:},'free_parameters',maxparams,'shocks',options_.ms.vd_shock_draws);  
-  if isscalar(vd)
-    mexErrCheck('mex_ms_variance_decomposition ergodic ',vd);
-  end
-  plot_ms_variance_decomposition(vd, 'Ergodic Variance Decomposition'); 
-  [regime_vd] = mex_ms_variance_decomposition(opt{:},'free_parameters',maxparams,'shocks',options_.ms.vd_shock_draws,'regimes');
-  if isscalar(regime_vd)
-    mexErrCheck('mex_ms_variance_decomposition ergodic regimes',regime_vd);
-  end
-  save([fname '/' fname '_ergodic_vd.mat'], 'vd', 'regime_vd');
-
-  if exist(mh_file,'file') > 0
-    [vd] = mex_ms_variance_decomposition(opt{:},'simulation_file',mh_file,'shocks',options_.ms.vd_shocks_per_parameter,'parameter_uncertainty');  
-    if isscalar(vd)
-      mexErrCheck('mex_ms_variance_decomposition bayesian ',vd);
-    end
-    [regime_vd] = mex_ms_variance_decomposition(opt{:},'simulation_file',mh_file,'shocks',options_.ms.vd_shocks_per_parameter,'parameter_uncertainty','regimes');
-    if isscalar(regime_vd)
-      mexErrCheck('mex_ms_variance_decomposition bayesian regimes ',regime_vd);
-    end
-    save([fname '/' fname '_bayesian_vd.mat'], 'vd', 'regime_vd');
-  end
-end
-
-%==========================================================================
diff --git a/ms-sbvar/Makefile b/ms-sbvar/Makefile
index 201a2b6917869e4087395363f6e18ada3a146e6a..475e85a0f71d9d19ef06a4fa8becb2411837175e 100644
--- a/ms-sbvar/Makefile
+++ b/ms-sbvar/Makefile
@@ -76,7 +76,7 @@ EXEC_OBJS = $(OBJS) dw_switch.o dw_switchio.o dw_dirichlet_restrictions.o dw_met
 # EXEC_OBJS := $(EXEC_OBJS)  tzmatlab.o csminwel.o
 
 # PROJECT FILES
-EXEC_OBJS := $(EXEC_OBJS) VARbase.o VARio.o dw_sbvar_command_line.o sbvar_estimate.o sbvar_simulate.o \
+EXEC_OBJS := $(EXEC_OBJS) VARbase.o VARio.o dw_sbvar_command_line.o sbvar_estimate.o sbvar_simulate.o sbvar_variance_decomposition.o \
                           sbvar_probabilities.o sbvar_mdd.o dw_csminwel.o sbvar_impulse_responses.o sbvar_forecast.o
 
 # OUTPUT
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index 5ab535a2674c1679968f4f22a9cc6a98f86352d4..c15009293b809cd35123e0f912f89fc7e8fcc29e 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -936,7 +936,7 @@ void
 MSSBVAREstimationStatement::writeOutput(ostream &output, const string &basename) const
 {
   options_list.writeOutput(output);
-  output << "ms_estimation(options_);" << endl;
+  output << "[options_, oo_] = ms_estimation(M_, options_, oo_);" << endl;
 }
 
 MSSBVARSimulationStatement::MSSBVARSimulationStatement(const OptionsList &options_list_arg) :
@@ -954,7 +954,7 @@ void
 MSSBVARSimulationStatement::writeOutput(ostream &output, const string &basename) const
 {
   options_list.writeOutput(output);
-  output << "ms_simulation(options_);" << endl;
+  output << "[options_, oo_] = ms_simulation(M_, options_, oo_);" << endl;
 }
 
 MSSBVARComputeMDDStatement::MSSBVARComputeMDDStatement(const OptionsList &options_list_arg) :
@@ -972,7 +972,7 @@ void
 MSSBVARComputeMDDStatement::writeOutput(ostream &output, const string &basename) const
 {
   options_list.writeOutput(output);
-  output << "oo_ = ms_compute_mdd(options_,oo_);" << endl;
+  output << "[options_, oo_] = ms_compute_mdd(M_, options_, oo_);" << endl;
 }
 
 MSSBVARComputeProbabilitiesStatement::MSSBVARComputeProbabilitiesStatement(const OptionsList &options_list_arg) :
@@ -990,7 +990,7 @@ void
 MSSBVARComputeProbabilitiesStatement::writeOutput(ostream &output, const string &basename) const
 {
   options_list.writeOutput(output);
-  output << "ms_compute_probabilities(options_);" << endl;
+  output << "[options_, oo_] = ms_compute_probabilities(M_, options_, oo_);" << endl;
 }
 
 MSSBVARIrfStatement::MSSBVARIrfStatement(const OptionsList &options_list_arg) :
@@ -1008,7 +1008,7 @@ void
 MSSBVARIrfStatement::writeOutput(ostream &output, const string &basename) const
 {
   options_list.writeOutput(output);
-  output << "ms_irf(options_);" << endl;
+  output << "[options_, oo_] = ms_irf(M_, options_, oo_);" << endl;
 }
 
 MSSBVARForecastStatement::MSSBVARForecastStatement(const OptionsList &options_list_arg) :
@@ -1026,7 +1026,7 @@ void
 MSSBVARForecastStatement::writeOutput(ostream &output, const string &basename) const
 {
   options_list.writeOutput(output);
-  output << "ms_forecast(options_);" << endl;
+  output << "[options_, oo_] = ms_forecast(M_, options_, oo_);" << endl;
 }
 
 MSSBVARVarianceDecompositionStatement::MSSBVARVarianceDecompositionStatement(const OptionsList &options_list_arg) :
@@ -1044,7 +1044,7 @@ void
 MSSBVARVarianceDecompositionStatement::writeOutput(ostream &output, const string &basename) const
 {
   options_list.writeOutput(output);
-  output << "ms_variance_decomposition(options_);" << endl;
+  output << "[options_, oo_] = ms_variance_decomposition(M_, options_, oo_);" << endl;
 }
 
 IdentificationStatement::IdentificationStatement(const OptionsList &options_list_arg)
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 736a242e732aa3ccf42cb3ec554210c5b398dbac..e6fe81d915003c4bf2471a2b78f9cea80003995f 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -151,16 +151,16 @@ class ParsingDriver;
 %token INDXPARR INDXOVR INDXAP APBAND INDXIMF IMFBAND INDXFORE FOREBAND INDXGFOREHAT INDXGIMFHAT
 %token INDXESTIMA INDXGDLS EQ_MS FILTER_COVARIANCE FILTER_DECOMPOSITION
 %token EQ_CMS TLINDX TLNUMBER BANACT CREATE_INITIALIZATION_FILE
-%token OUTPUT_FILE_TAG DRAWS_NBR_BURN_IN_1 DRAWS_NBR_BURN_IN_2
+%token OUTPUT_FILE_TAG DRAWS_NBR_BURN_IN_1 DRAWS_NBR_BURN_IN_2 HORIZON
 %token SBVAR TREND_VAR DEFLATOR GROWTH_FACTOR MS_IRF MS_VARIANCE_DECOMPOSITION
 %token MS_ESTIMATION MS_SIMULATION MS_COMPUTE_MDD MS_COMPUTE_PROBABILITIES MS_FORECAST
 %token SVAR_IDENTIFICATION EQUATION EXCLUSION LAG UPPER_CHOLESKY LOWER_CHOLESKY
 %token MARKOV_SWITCHING CHAIN STATE DURATION NUMBER_OF_STATES
 %token SVAR COEFFICIENTS VARIANCES CONSTANTS EQUATIONS
 %token EXTERNAL_FUNCTION EXT_FUNC_NAME EXT_FUNC_NARGS FIRST_DERIV_PROVIDED SECOND_DERIV_PROVIDED
-%token SELECTED_VARIABLES_ONLY COVA_COMPUTE
-%token ERROR_BANDS ERROR_BAND_PERCENTILES PARAMETER_UNCERTAINTY
-%token SHOCK_DRAWS REGIMES FREE_PARAMETERS MEDIAN DATA_OBS_NBR
+%token SELECTED_VARIABLES_ONLY COVA_COMPUTE ESTIMATION_FILE_TAG SIMULATION_FILE_TAG
+%token ERROR_BANDS ERROR_BAND_PERCENTILES SHOCKS_PER_PARAMETER
+%token SHOCK_DRAWS FREE_PARAMETERS MEDIAN DATA_OBS_NBR
 %token FILTERED_PROBABILITIES FILTERED REAL_TIME_SMOOTHED
 %token PROPOSAL_TYPE MDD_PROPOSAL_DRAWS MDD_USE_MEAN_CENTER
 %token ADAPTIVE_MH_DRAWS THINNING_FACTOR COEFFICIENTS_PRIOR_HYPERPARAMETERS
@@ -1473,13 +1473,14 @@ sbvar : SBVAR ';'
       ;
 
 ms_variance_decomposition_option : o_output_file_tag
+                                 | o_estimation_file_tag
+                                 | o_simulation_file_tag
                                  | o_filtered_probabilities
                                  | o_error_bands
                                  | o_error_band_percentiles
-                                 | o_parameter_uncertainty
                                  | o_shock_draws
+                                 | o_shocks_per_parameter
                                  | o_thinning_factor
-                                 | o_regimes
                                  | o_free_parameters
                                  | o_load_mh_file
                                  | o_median
@@ -1496,13 +1497,14 @@ ms_variance_decomposition : MS_VARIANCE_DECOMPOSITION ';'
                           ;
 
 ms_forecast_option : o_output_file_tag
+                   | o_estimation_file_tag
+                   | o_simulation_file_tag
                    | o_data_obs_nbr
                    | o_error_bands
                    | o_error_band_percentiles
-                   | o_parameter_uncertainty
                    | o_shock_draws
+                   | o_shocks_per_parameter
                    | o_thinning_factor
-                   | o_regimes
                    | o_free_parameters
                    | o_load_mh_file
                    | o_median
@@ -1519,13 +1521,15 @@ ms_forecast : MS_FORECAST ';'
             ;
 
 ms_irf_option : o_output_file_tag
+              | o_estimation_file_tag
+              | o_simulation_file_tag
+              | o_horizon
               | o_filtered_probabilities
               | o_error_bands
               | o_error_band_percentiles
-              | o_parameter_uncertainty
               | o_shock_draws
+              | o_shocks_per_parameter
               | o_thinning_factor
-              | o_regimes
               | o_free_parameters
               | o_load_mh_file
               | o_median
@@ -1542,6 +1546,8 @@ ms_irf : MS_IRF ';'
        ;
 
 ms_compute_probabilities_option : o_output_file_tag
+                                | o_estimation_file_tag
+                                | o_simulation_file_tag
                                 | o_filtered_probabilities
                                 | o_real_time_smoothed
                                 ;
@@ -1557,6 +1563,8 @@ ms_compute_probabilities : MS_COMPUTE_PROBABILITIES ';'
                          ;
 
 ms_compute_mdd_option : o_output_file_tag
+                      | o_estimation_file_tag
+                      | o_simulation_file_tag
                       | o_load_mh_file
                       | o_proposal_type
                       | o_mdd_proposal_draws
@@ -1574,6 +1582,7 @@ ms_compute_mdd : MS_COMPUTE_MDD ';'
                ;
 
 ms_simulation_option : o_output_file_tag
+                     | o_estimation_file_tag
                      | o_mh_replic
                      | o_drop
                      | o_thinning_factor
@@ -1947,10 +1956,8 @@ o_parameter_set : PARAMETER_SET EQUAL PRIOR_MODE
                 | PARAMETER_SET EQUAL POSTERIOR_MEDIAN
                   { driver.option_str("parameter_set", "posterior_median"); }
                 ;
-
 o_shocks : SHOCKS EQUAL '(' list_of_symbol_lists ')' { driver.option_symbol_list("shocks"); };
 o_labels : LABELS EQUAL '(' symbol_list ')' { driver.option_symbol_list("labels"); };
-
 o_freq : FREQ EQUAL INT_NUMBER {driver.option_num("ms.freq",$3); };
 o_initial_year : INITIAL_YEAR EQUAL INT_NUMBER {driver.option_num("ms.initial_year",$3); };
 o_initial_subperiod : INITIAL_SUBPERIOD EQUAL INT_NUMBER {driver.option_num("ms.initial_subperiod",$3); };
@@ -2044,6 +2051,8 @@ o_selected_variables_only : SELECTED_VARIABLES_ONLY
 o_cova_compute : COVA_COMPUTE EQUAL INT_NUMBER
                  { driver.option_num("cova_compute",$3);}
                ;
+o_estimation_file_tag : ESTIMATION_FILE_TAG EQUAL filename { driver.option_str("ms.estimation_file_tag", $3); };
+o_simulation_file_tag : SIMULATION_FILE_TAG EQUAL filename { driver.option_str("ms.simulation_file_tag", $3); };
 o_upper_cholesky : UPPER_CHOLESKY { driver.option_num("ms.upper_cholesky","1"); };
 o_lower_cholesky : LOWER_CHOLESKY { driver.option_num("ms.lower_cholesky","1"); };
 o_coefficients_prior_hyperparameters : COEFFICIENTS_PRIOR_HYPERPARAMETERS EQUAL vec_value
@@ -2083,14 +2092,13 @@ o_adaptive_mh_draws : ADAPTIVE_MH_DRAWS EQUAL INT_NUMBER { driver.option_num("ms
 o_mdd_proposal_draws : MDD_PROPOSAL_DRAWS EQUAL INT_NUMBER { driver.option_num("ms.mdd_proposal_draws",$3); };
 o_mdd_use_mean_center : MDD_USE_MEAN_CENTER { driver.option_num("ms.mdd_use_mean_center","1"); };
 o_proposal_type : PROPOSAL_TYPE EQUAL vec_value { driver.option_num("ms.proposal_type",$3); };
+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"); };
-
 o_error_bands : ERROR_BANDS { driver.option_num("ms.error_bands","1"); };
 o_error_band_percentiles : ERROR_BAND_PERCENTILES EQUAL vec_value { driver.option_num("ms.percentiles",$3); };
-o_parameter_uncertainty : PARAMETER_UNCERTAINTY { driver.option_num("ms.error_bands","1"); };
 o_shock_draws : SHOCK_DRAWS EQUAL INT_NUMBER { driver.option_num("ms.shock_draws",$3); };
-o_regimes : REGIMES { driver.option_num("ms.regimes","1"); };
+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_data_obs_nbr : DATA_OBS_NBR { driver.option_num("ms.forecast_data_obs","1"); };
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index d1a4b452de921ecd198275fa84258294d4bbba8d..dd988128a9d190494ee047e83aee604c871e5aed 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -320,13 +320,15 @@ string eofbuff;
 <DYNARE_STATEMENT>banact {return token::BANACT;}
 
 <DYNARE_STATEMENT>output_file_tag {return token::OUTPUT_FILE_TAG;}
+<DYNARE_STATEMENT>estimation_file_tag {return token::ESTIMATION_FILE_TAG;};
+<DYNARE_STATEMENT>simulation_file_tag {return token::SIMULATION_FILE_TAG;};
 <DYNARE_STATEMENT>filtered {return token::FILTERED;}
+<DYNARE_STATEMENT>horizon {return token::HORIZON;}
 <DYNARE_STATEMENT>error_bands {return token::ERROR_BANDS;}
 <DYNARE_STATEMENT>error_band_percentiles {return token::ERROR_BAND_PERCENTILES;}
-<DYNARE_STATEMENT>parameter_uncertainty {return token::PARAMETER_UNCERTAINTY;}
 <DYNARE_STATEMENT>shock_draws {return token::SHOCK_DRAWS;}
+<DYNARE_STATEMENT>shocks_per_parameter {return token::SHOCKS_PER_PARAMETER;}
 <DYNARE_STATEMENT>thinning_factor {return token::THINNING_FACTOR;}
-<DYNARE_STATEMENT>regimes {return token::REGIMES;}
 <DYNARE_STATEMENT>free_parameters {return token::FREE_PARAMETERS;}
 <DYNARE_STATEMENT>median {return token::MEDIAN;}
 <DYNARE_STATEMENT>data_obs_nbr {return token::DATA_OBS_NBR;}
diff --git a/tests/ms-sbvar/test_ms_variances.mod b/tests/ms-sbvar/test_ms_variances.mod
index 1380683def1f7adc1410ac44ef810686b148007a..65d1997a3455b9c46a29c8bf49d1dbe528aac9fb 100644
--- a/tests/ms-sbvar/test_ms_variances.mod
+++ b/tests/ms-sbvar/test_ms_variances.mod
@@ -16,7 +16,12 @@ markov_switching(chain=1,number_of_states=2,duration=2.5);
 
 svar(variances, chain=1);
 
+set_dynare_seed(5);
+
 ms_estimation(datafile=data,freq=4,initial_year=1959,final_year=2005,nlags=4);
 ms_simulation(mh_replic=1000);
 ms_compute_mdd;
 ms_compute_probabilities;
+ms_irf;
+ms_forecast;
+ms_variance_decomposition;