diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index ed6abf87d96f0b1a9de4cc1bc2538e6245773461..1cfb81abd9655f35836cfc555db206dab48d525b 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -255,14 +255,13 @@ options_.initial_date.subperiod = 0;
 % discretionary policy
 options_.discretionary_policy = 0;
 
-%% SWZ SBVAR
-options_.ms.filtered_probabilities = 1;
-options_.ms.real_time_smoothed_probabilities = 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.restriction_fname = 0;
 options_.ms.cross_restrictions = 0;
 options_.ms.contemp_reduced_form = 0;
 options_.ms.bayesian_prior = 1;
@@ -272,38 +271,21 @@ options_.ms.gsig2_lmd = 50^2;
 options_.ms.gsig2_lmdm = 50^2;
 options_.ms.lower_cholesky = 0;
 options_.ms.upper_cholesky = 0;
-options_.ms.create_initialization_file = 1;
-options_.ms.compute_mdd = 1;
-options_.ms.compute_probabilities = 1;
-options_.ms.coefficients_prior_hyperparameters = [1.0 1.0 0.1 1.2 1.0 1.0];
-options_.ms.irf = 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.forecast = 0;
-options_.ms.variance_decomposition = 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.estimate.random = 5;
-options_.ms.estimate.random_max = 20;
-options_.ms.estimate.random_tol_obj = 0.1;
-options_.ms.estimate.random_tol_parms = 0.1;
-options_.ms.estimate.cb = 1e-3;
-options_.ms.estimate.ce = 1e-6;
-options_.ms.estimate.ci = 0.1;
-options_.ms.estimate.ib = 50;
-options_.ms.estimate.ii = 2.0;
-options_.ms.estimate.mb = 100;
 options_.ms.load_mh_file = 0;
-options_.ms.mh_replic = 20000;
-options_.ms.thinning_factor = 10;
-options_.ms.drop = 2000;
-options_.ms.adapative_mh_draws = 30000;
-options_.ms.mdd_proposal_type = [3 0.1 0.9];
-options_.ms.mdd_use_mean_center = 0;
-options_.ms.mdd_proposal_draws = 100000;
 options_.ms.filtered_probabilities = 0;
 options_.ms.real_time_smoothed_probabilities = 0;
 options_.ms.irf_shocks_per_parameter = options_.ms.shocks_per_parameter;
diff --git a/matlab/ms-sbvar/clean_ms_files.m b/matlab/ms-sbvar/clean_ms_files.m
index 2d0e86d3544a3f9c6fef54b0cc3d5d95767b2135..a19fee68e81d60523879e390d9ed50f8528c9e2e 100644
--- a/matlab/ms-sbvar/clean_ms_files.m
+++ b/matlab/ms-sbvar/clean_ms_files.m
@@ -1,4 +1,4 @@
-function clean_ms_files(mod_name)
+function clean_ms_files(output_file_tag)
 % function clean_ms_files()
 % removes MS intermediary files
 %
@@ -28,29 +28,26 @@ function clean_ms_files(mod_name)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-    delete_if_exist(['./draws_' mod_name '.dat'])
-    delete_if_exist(['est_aux_' mod_name '.dat'])
-    delete_if_exist(['est_csminwel_' mod_name '.dat'])
-    delete_if_exist(['est_final_' mod_name '.dat'])
-    delete_if_exist(['est_flat_header_' mod_name '.dat'])
-    delete_if_exist(['est_flat_' mod_name '.dat'])
-    delete_if_exist(['est_intermediate_' mod_name '.dat'])
-    delete_if_exist(['header_' mod_name '.dat'])
-    delete_if_exist(['init_' mod_name '.dat'])
-    delete_if_exist('markov_file.dat')
-    delete_if_exist(['matlab_' mod_name '.prn'])
-    delete_if_exist(['mhm_draws_' mod_name '.dat'])
-    delete_if_exist(['mhm_input_' mod_name '.dat'])
-    delete_if_exist(['mhm_intermediate_draws_' mod_name '.dat'])
-    delete_if_exist(['mhm_intermediate_' mod_name '.dat'])
-    delete_if_exist(['mhm_regime_counts_' mod_name '.dat'])
-    delete_if_exist(['probabilities_' mod_name '.dat'])
-    delete_if_exist(['truncatedpower_md_posterior_' mod_name '.dat'])
-    delete_if_exist(['truncatedpower_md_proposal_' mod_name '.dat'])
-    delete_if_exist(['truncatedpower_md_' mod_name '.dat'])
+delete_if_exist(['est_aux_' output_file_tag '.out']);
+delete_if_exist(['est_csminwel_' output_file_tag '.out']);
+delete_if_exist(['est_final_' output_file_tag '.out']);
+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('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_t3_' output_file_tag '.out']);
+delete_if_exist(['proposal_t3_' output_file_tag '.out']);
+delete_if_exist(['simulation_info_' output_file_tag '.out']);
+delete_if_exist(['simulation_' output_file_tag '.out']);
+delete_if_exist([output_file_tag '_markov_file.dat']);
+end
 
 function delete_if_exist(fname)
-    if exist(fname) == 2
-        delete(fname)
-    end
-    
\ No newline at end of file
+if exist(fname,'file') == 2
+    delete(fname);
+end
+end
diff --git a/matlab/ms-sbvar/ms_compute_mdd.m b/matlab/ms-sbvar/ms_compute_mdd.m
new file mode 100644
index 0000000000000000000000000000000000000000..387c89737103c3dff9608ae5721e1d8ecb4b4e22
--- /dev/null
+++ b/matlab/ms-sbvar/ms_compute_mdd.m
@@ -0,0 +1,78 @@
+function oo_=ms_compute_mdd(options_,oo_)
+%function ms_compute_mdd()
+% calls ms sbvar mdd mex file
+%
+% INPUTS
+%    options_
+%
+% 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/>.
+
+compute_mdd = create_mdd_commandline(options_);
+[err] = ms_sbvar_command_line(compute_mdd);
+mexErrCheck('ms_sbvar_command_line mdd',err);
+
+% 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
+
+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
+end
diff --git a/matlab/ms-sbvar/ms_compute_probabilities.m b/matlab/ms-sbvar/ms_compute_probabilities.m
new file mode 100644
index 0000000000000000000000000000000000000000..42b6b07d8c88b72a73c1d49982a571e3c5473429
--- /dev/null
+++ b/matlab/ms-sbvar/ms_compute_probabilities.m
@@ -0,0 +1,47 @@
+function ms_compute_probabilities(options_)
+%function ms_simulation()
+% calls ms sbvar probabilities mex file
+%
+% INPUTS
+%    options_
+%
+% 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/>.
+
+compute_probabilities = create_probabilities_commandline(options_);
+[err] = ms_sbvar_command_line(compute_probabilities);
+mexErrCheck('ms_sbvar_command_line probabilities',err);
+
+end
+
+function opt=create_probabilities_commandline(options_)
+
+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')];
+
+if options_.DynareRandomStreams.seed
+    opt = [' -seed' num2str(options_.DynareRandomStreams.seed)];
+end
+end
diff --git a/matlab/ms-sbvar/ms_estimation.m b/matlab/ms-sbvar/ms_estimation.m
new file mode 100644
index 0000000000000000000000000000000000000000..400e62773aebe4175e0fcf5787c5cce2cfffa365
--- /dev/null
+++ b/matlab/ms-sbvar/ms_estimation.m
@@ -0,0 +1,72 @@
+function ms_estimation(options_)
+%function ms_estimation()
+% calls ms sbvar estimation mex file
+%
+% INPUTS
+%    options_
+%
+% 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/>.
+
+% 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
+
+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
+
+end
diff --git a/matlab/ms-sbvar/ms_sbvar_setup.m b/matlab/ms-sbvar/ms_sbvar_setup.m
new file mode 100644
index 0000000000000000000000000000000000000000..cefcdd93556214fa152c453e9b4e1c2932f27f6b
--- /dev/null
+++ b/matlab/ms-sbvar/ms_sbvar_setup.m
@@ -0,0 +1,453 @@
+function ms_sbvar_setup(options_)
+%function ms_sbvar_setup()
+% does the general file initialization for ms sbvar
+%
+% INPUTS
+%    options_
+%
+% OUTPUTS
+%    none
+%
+% SPECIAL REQUIREMENTS
+%    none
+
+% Copyright (C) 2003-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/>.
+
+dynareroot = strrep(which('dynare'),'dynare.m','');
+ms_root = [dynareroot '/ms-sbvar'];
+addpath([ms_root '/cstz']);
+addpath([ms_root '/identification']);
+addpath([ms_root '/switching_specification']);
+addpath([ms_root '/mhm_specification']);
+
+if ~isfield(options_.ms,'initial_year')
+    error('Must set initial_year option');
+end
+
+if ~isfield(options_.ms,'final_year')
+    error('Must set final_year option');
+end
+
+if ~isfield(options_,'datafile')
+    error('Must set datafile option');
+end
+
+options_.data = read_variables(options_.datafile, ...
+    options_.varobs, [], options_.xls_sheet, options_.xls_range);
+
+if ~isfield(options_.ms, 'forecast')
+    options_.ms.forecast = 4;
+end
+
+if options_.ms.upper_cholesky
+    if options_.ms.lower_cholesky
+        error(['Upper Cholesky and lower Cholesky decomposition can''t be ' ...
+               'requested at the same time!'])
+    else
+        options_.ms.restriction_fname = 'upper_cholesky';
+    end
+elseif options_.ms.lower_cholesky
+    options_.ms.restriction_fname = 'lower_cholesky';
+%elseif ~isempty(options_.ms.Qi) && ~isempty(options_.ms.Ri)
+%    options_.ms.restriction_fname = 'exclusions';
+else
+    options_.ms.restriction_fname = 0;
+end
+
+%==========================================================================
+%== Markov Process Specification File
+%==========================================================================
+markov_file = [options_.ms.output_file_tag '_markov_file.dat'];
+
+%==========================================================================
+%== 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 = options_.ms.alpha;
+
+% Beta on p. 66 for squared time-varying structural shock lambda.
+gbeta = options_.ms.beta;
+
+% 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 = options_.ms.gsig2_lmdm;
+
+
+%==========================================================================
+%== Data
+%==========================================================================
+% Read in data to produce rectangular array named xdd.  Each column is one
+% data series.
+xdd=options_.data;
+
+% Information about timing of the data for consistancy checks
+% quarters (4) or months (12)
+q_m = options_.ms.freq;
+% beginning year in data set
+yrBin=options_.ms.initial_year;
+% beginning quarter or month in data set
+%options_.ms.initial_subperiod = 1;
+qmBin=options_.ms.initial_subperiod;
+% final year in data set
+yrFin=options_.ms.final_year;
+% final month or quarter in data set
+qmFin=options_.ms.final_subperiod;
+% first year to use in estimation
+yrStart=options_.ms.initial_year;
+% first quarter or month to use in estimation
+qmStart=options_.ms.initial_subperiod;
+% last year to use in estimation
+yrEnd=options_.ms.final_year;
+% last quater or month to use in estimation
+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;
+
+%==========================================================================
+%==========================================================================
+%==========================================================================
+%== 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
+%==========================================================================
+
+%======================================================================
+%== 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,options_)
+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
\ No newline at end of file
diff --git a/matlab/ms-sbvar/ms_simulation.m b/matlab/ms-sbvar/ms_simulation.m
new file mode 100644
index 0000000000000000000000000000000000000000..adbc0cb5812cee0187cf581a3b70980647519b49
--- /dev/null
+++ b/matlab/ms-sbvar/ms_simulation.m
@@ -0,0 +1,49 @@
+function ms_simulation(options_)
+%function ms_simulation()
+% calls ms sbvar simulation mex file
+%
+% INPUTS
+%    options_
+%
+% 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/>.
+
+perform_simulation = create_simulation_commandline(options_);
+[err] = ms_sbvar_command_line(perform_simulation);
+mexErrCheck('ms_sbvar_command_line simulate',err);
+
+end
+
+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')];
+
+if options_.DynareRandomStreams.seed
+    opt = [' -seed' num2str(options_.DynareRandomStreams.seed)];
+end
+end
diff --git a/matlab/ms-sbvar/ms_write_markov_file.m b/matlab/ms-sbvar/ms_write_markov_file.m
index 2f04bee5b73708d3f6b911bce82914c26c592ab0..b2bcc5d43caf765b01f0156cc62783268af4675c 100644
--- a/matlab/ms-sbvar/ms_write_markov_file.m
+++ b/matlab/ms-sbvar/ms_write_markov_file.m
@@ -1,4 +1,4 @@
-function ms_write_markov_file(fname,M,options)
+function ms_write_markov_file(fname,options)
 
     n_chains = length(options.ms.ms_chain);
     nvars = size(options.varobs,1);
diff --git a/matlab/ms-sbvar/set_flag_if_exists.m b/matlab/ms-sbvar/set_flag_if_exists.m
new file mode 100644
index 0000000000000000000000000000000000000000..cf6f2e24fba39791c9b9cb437b4b1c6e919e1bca
--- /dev/null
+++ b/matlab/ms-sbvar/set_flag_if_exists.m
@@ -0,0 +1,97 @@
+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/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index ed1454ed7ee52c5d5ced8dae7eedde566ff5c22f..5ab535a2674c1679968f4f22a9cc6a98f86352d4 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -921,22 +921,130 @@ SBVARStatement::writeOutput(ostream &output, const string &basename) const
   output << "ms_sbvar(0,M_,options_);" << endl;
 }
 
-MS_SBVARStatement::MS_SBVARStatement(const OptionsList &options_list_arg) :
+MSSBVAREstimationStatement::MSSBVAREstimationStatement(const OptionsList &options_list_arg) :
   options_list(options_list_arg)
 {
 }
 
 void
-MS_SBVARStatement::checkPass(ModFileStructure &mod_file_struct)
+MSSBVAREstimationStatement::checkPass(ModFileStructure &mod_file_struct)
 {
   mod_file_struct.bvar_present = true;
 }
 
 void
-MS_SBVARStatement::writeOutput(ostream &output, const string &basename) const
+MSSBVAREstimationStatement::writeOutput(ostream &output, const string &basename) const
 {
   options_list.writeOutput(output);
-  output << "ms_sbvar(1,M_,options_);" << endl;
+  output << "ms_estimation(options_);" << endl;
+}
+
+MSSBVARSimulationStatement::MSSBVARSimulationStatement(const OptionsList &options_list_arg) :
+  options_list(options_list_arg)
+{
+}
+
+void
+MSSBVARSimulationStatement::checkPass(ModFileStructure &mod_file_struct)
+{
+  mod_file_struct.bvar_present = true;
+}
+
+void
+MSSBVARSimulationStatement::writeOutput(ostream &output, const string &basename) const
+{
+  options_list.writeOutput(output);
+  output << "ms_simulation(options_);" << endl;
+}
+
+MSSBVARComputeMDDStatement::MSSBVARComputeMDDStatement(const OptionsList &options_list_arg) :
+  options_list(options_list_arg)
+{
+}
+
+void
+MSSBVARComputeMDDStatement::checkPass(ModFileStructure &mod_file_struct)
+{
+  mod_file_struct.bvar_present = true;
+}
+
+void
+MSSBVARComputeMDDStatement::writeOutput(ostream &output, const string &basename) const
+{
+  options_list.writeOutput(output);
+  output << "oo_ = ms_compute_mdd(options_,oo_);" << endl;
+}
+
+MSSBVARComputeProbabilitiesStatement::MSSBVARComputeProbabilitiesStatement(const OptionsList &options_list_arg) :
+  options_list(options_list_arg)
+{
+}
+
+void
+MSSBVARComputeProbabilitiesStatement::checkPass(ModFileStructure &mod_file_struct)
+{
+  mod_file_struct.bvar_present = true;
+}
+
+void
+MSSBVARComputeProbabilitiesStatement::writeOutput(ostream &output, const string &basename) const
+{
+  options_list.writeOutput(output);
+  output << "ms_compute_probabilities(options_);" << endl;
+}
+
+MSSBVARIrfStatement::MSSBVARIrfStatement(const OptionsList &options_list_arg) :
+  options_list(options_list_arg)
+{
+}
+
+void
+MSSBVARIrfStatement::checkPass(ModFileStructure &mod_file_struct)
+{
+  mod_file_struct.bvar_present = true;
+}
+
+void
+MSSBVARIrfStatement::writeOutput(ostream &output, const string &basename) const
+{
+  options_list.writeOutput(output);
+  output << "ms_irf(options_);" << endl;
+}
+
+MSSBVARForecastStatement::MSSBVARForecastStatement(const OptionsList &options_list_arg) :
+  options_list(options_list_arg)
+{
+}
+
+void
+MSSBVARForecastStatement::checkPass(ModFileStructure &mod_file_struct)
+{
+  mod_file_struct.bvar_present = true;
+}
+
+void
+MSSBVARForecastStatement::writeOutput(ostream &output, const string &basename) const
+{
+  options_list.writeOutput(output);
+  output << "ms_forecast(options_);" << endl;
+}
+
+MSSBVARVarianceDecompositionStatement::MSSBVARVarianceDecompositionStatement(const OptionsList &options_list_arg) :
+  options_list(options_list_arg)
+{
+}
+
+void
+MSSBVARVarianceDecompositionStatement::checkPass(ModFileStructure &mod_file_struct)
+{
+  mod_file_struct.bvar_present = true;
+}
+
+void
+MSSBVARVarianceDecompositionStatement::writeOutput(ostream &output, const string &basename) const
+{
+  options_list.writeOutput(output);
+  output << "ms_variance_decomposition(options_);" << endl;
 }
 
 IdentificationStatement::IdentificationStatement(const OptionsList &options_list_arg)
diff --git a/preprocessor/ComputingTasks.hh b/preprocessor/ComputingTasks.hh
index 4c22c14004d70752dea2c7298ba051a0564660bf..c05fa545e991777c15dbcd02e509ff829d40a331 100644
--- a/preprocessor/ComputingTasks.hh
+++ b/preprocessor/ComputingTasks.hh
@@ -378,12 +378,72 @@ public:
   virtual void writeOutput(ostream &output, const string &basename) const;
 };
 
-class MS_SBVARStatement : public Statement
+class MSSBVAREstimationStatement : public Statement
 {
 private:
   const OptionsList options_list;
 public:
-  MS_SBVARStatement(const OptionsList &options_list_arg);
+  MSSBVAREstimationStatement(const OptionsList &options_list_arg);
+  virtual void checkPass(ModFileStructure &mod_file_struct);
+  virtual void writeOutput(ostream &output, const string &basename) const;
+};
+
+class MSSBVARSimulationStatement : public Statement
+{
+private:
+  const OptionsList options_list;
+public:
+  MSSBVARSimulationStatement(const OptionsList &options_list_arg);
+  virtual void checkPass(ModFileStructure &mod_file_struct);
+  virtual void writeOutput(ostream &output, const string &basename) const;
+};
+
+class MSSBVARComputeMDDStatement : public Statement
+{
+private:
+  const OptionsList options_list;
+public:
+  MSSBVARComputeMDDStatement(const OptionsList &options_list_arg);
+  virtual void checkPass(ModFileStructure &mod_file_struct);
+  virtual void writeOutput(ostream &output, const string &basename) const;
+};
+
+class MSSBVARComputeProbabilitiesStatement : public Statement
+{
+private:
+  const OptionsList options_list;
+public:
+  MSSBVARComputeProbabilitiesStatement(const OptionsList &options_list_arg);
+  virtual void checkPass(ModFileStructure &mod_file_struct);
+  virtual void writeOutput(ostream &output, const string &basename) const;
+};
+
+class MSSBVARIrfStatement : public Statement
+{
+private:
+  const OptionsList options_list;
+public:
+  MSSBVARIrfStatement(const OptionsList &options_list_arg);
+  virtual void checkPass(ModFileStructure &mod_file_struct);
+  virtual void writeOutput(ostream &output, const string &basename) const;
+};
+
+class MSSBVARForecastStatement : public Statement
+{
+private:
+  const OptionsList options_list;
+public:
+  MSSBVARForecastStatement(const OptionsList &options_list_arg);
+  virtual void checkPass(ModFileStructure &mod_file_struct);
+  virtual void writeOutput(ostream &output, const string &basename) const;
+};
+
+class MSSBVARVarianceDecompositionStatement : public Statement
+{
+private:
+  const OptionsList options_list;
+public:
+  MSSBVARVarianceDecompositionStatement(const OptionsList &options_list_arg);
   virtual void checkPass(ModFileStructure &mod_file_struct);
   virtual void writeOutput(ostream &output, const string &basename) const;
 };
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index ee77a818d2e2b5f8520fa685db33b4ddeda655aa..736a242e732aa3ccf42cb3ec554210c5b398dbac 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -150,16 +150,26 @@ class ParsingDriver;
 %token GSIG2_LMD GSIG2_LMDM Q_DIAG FLAT_PRIOR NCSK NSTD
 %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 ESTIMATE_MSMODEL
-%token COMPUTE_MDD COMPUTE_PROBABILITIES PRINT_DRAWS N_DRAWS THINNING_FACTOR PROPOSAL_DRAWS MARKOV_FILE
-%token MHM_FILE OUTPUT_FILE_TAG DRAWS_NBR_BURN_IN_1 DRAWS_NBR_BURN_IN_2 DRAWS_NBR_MEAN_VAR_ESTIMATE
-%token DRAWS_NBR_MODIFIED_HARMONIC_MEAN DIRICHLET_SCALE
-%token SBVAR MS_SBVAR TREND_VAR DEFLATOR GROWTH_FACTOR
+%token EQ_CMS TLINDX TLNUMBER BANACT CREATE_INITIALIZATION_FILE
+%token OUTPUT_FILE_TAG DRAWS_NBR_BURN_IN_1 DRAWS_NBR_BURN_IN_2
+%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 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
+%token CONVERGENCE_STARTING_VALUE CONVERGENCE_ENDING_VALUE CONVERGENCE_INCREMENT_VALUE
+%token MAX_ITERATIONS_STARTING_VALUE MAX_ITERATIONS_INCREMENT_VALUE MAX_BLOCK_ITERATIONS
+%token MAX_REPEATED_OPTIMIZATION_RUNS FUNCTION_CONVERGENCE_CRITERION
+%token PARAMETER_CONVERGENCE_CRITERION NUMBER_OF_LARGE_PERTURBATIONS NUMBER_OF_SMALL_PERTURBATIONS
+%token NUMBER_OF_POSTERIOR_DRAWS_AFTER_PERTURBATION MAX_NUMBER_OF_STAGES
+%token RANDOM_FUNCTION_CONVERGENCE_CRITERION RANDOM_PARAMETER_CONVERGENCE_CRITERION
 
 %type <node_val> expression expression_or_empty
 %type <node_val> equation hand_side
@@ -221,7 +231,6 @@ statement : parameters
           | bvar_density
           | bvar_forecast
           | sbvar
-          | ms_sbvar
           | dynare_sensitivity
           | homotopy_setup
           | forecast
@@ -240,6 +249,13 @@ statement : parameters
           | external_function
           | steady_state_model
           | trend_var
+          | ms_estimation
+          | ms_simulation
+          | ms_compute_mdd
+          | ms_compute_probabilities
+          | ms_forecast
+          | ms_irf
+          | ms_variance_decomposition
           ;
 
 dsample : DSAMPLE INT_NUMBER ';'
@@ -1456,83 +1472,169 @@ sbvar : SBVAR ';'
         { driver.sbvar(); }
       ;
 
-ms_sbvar_option : o_datafile
-                | o_freq
-                | o_initial_year
-                | o_initial_subperiod
-                | o_final_year
-                | o_final_subperiod
-                | o_data
-                | o_vlist
-                | o_vlistlog
-                | o_vlistper
-                | o_varlist
-                | o_restriction_fname
-                | o_nlags
-                | o_cross_restrictions
-                | o_contemp_reduced_form
-                | o_real_pseudo_forecast
-                | o_bayesian_prior
-                | o_dummy_obs
-                | o_nstates
-                | o_indxscalesstates
-                | o_alpha
-                | o_beta
-                | o_gsig2_lmd
-                | o_gsig2_lmdm
-                | o_q_diag
-                | o_flat_prior
-                | o_ncsk
-                | o_nstd
-                | o_ninv
-                | o_indxparr
-                | o_indxovr
-                | o_aband
-                | o_indxap
-                | o_apband
-                | o_indximf
-                | o_indxfore
-                | o_foreband
-                | o_indxgforhat
-                | o_indxgimfhat
-                | o_indxestima
-                | o_indxgdls
-                | o_eq_ms
-                | o_cms
-                | o_ncms
-                | o_eq_cms
-                | o_tlindx
-                | o_tlnumber
-                | o_cnum
-                | o_forecast
-                | o_output_file_tag
-                | o_create_initialization_file
-                | o_estimate_msmodel
-                | o_compute_mdd
-                | o_compute_probabilities
-                | o_print_draws
-                | o_n_draws
-                | o_thinning_factor
-                | o_markov_file
-                | o_mhm_file
-                | o_proposal_draws
-                | o_draws_nbr_burn_in_1
-                | o_draws_nbr_burn_in_2
-                | o_draws_nbr_mean_var_estimate
-                | o_draws_nbr_modified_harmonic_mean
-                | o_dirichlet_scale
-                ;
+ms_variance_decomposition_option : o_output_file_tag
+                                 | o_filtered_probabilities
+                                 | o_error_bands
+                                 | o_error_band_percentiles
+                                 | o_parameter_uncertainty
+                                 | o_shock_draws
+                                 | o_thinning_factor
+                                 | o_regimes
+                                 | o_free_parameters
+                                 | o_load_mh_file
+                                 | o_median
+                                 ;
+
+ms_variance_decomposition_options_list : ms_variance_decomposition_option COMMA ms_variance_decomposition_options_list
+                                       | ms_variance_decomposition_option
+                                       ;
+
+ms_variance_decomposition : MS_VARIANCE_DECOMPOSITION ';'
+                            { driver.ms_variance_decomposition(); }
+                          | MS_VARIANCE_DECOMPOSITION '(' ms_variance_decomposition_options_list ')' ';'
+                            { driver.ms_variance_decomposition(); }
+                          ;
+
+ms_forecast_option : o_output_file_tag
+                   | o_data_obs_nbr
+                   | o_error_bands
+                   | o_error_band_percentiles
+                   | o_parameter_uncertainty
+                   | o_shock_draws
+                   | o_thinning_factor
+                   | o_regimes
+                   | o_free_parameters
+                   | o_load_mh_file
+                   | o_median
+                   ;
+
+ms_forecast_options_list : ms_forecast_option COMMA ms_forecast_options_list
+                         | ms_forecast_option
+                         ;
+
+ms_forecast : MS_FORECAST ';'
+              { driver.ms_forecast(); }
+            | MS_FORECAST '(' ms_forecast_options_list ')' ';'
+              { driver.ms_forecast(); }
+            ;
+
+ms_irf_option : o_output_file_tag
+              | o_filtered_probabilities
+              | o_error_bands
+              | o_error_band_percentiles
+              | o_parameter_uncertainty
+              | o_shock_draws
+              | o_thinning_factor
+              | o_regimes
+              | o_free_parameters
+              | o_load_mh_file
+              | o_median
+              ;
+
+ms_irf_options_list : ms_irf_option COMMA ms_irf_options_list
+                    | ms_irf_option
+                    ;
+
+ms_irf : MS_IRF ';'
+         { driver.ms_irf(); }
+       | MS_IRF '(' ms_irf_options_list ')' ';'
+         { driver.ms_irf(); }
+       ;
+
+ms_compute_probabilities_option : o_output_file_tag
+                                | o_filtered_probabilities
+                                | o_real_time_smoothed
+                                ;
 
-ms_sbvar_options_list : ms_sbvar_option COMMA ms_sbvar_options_list
-                      | ms_sbvar_option
+ms_compute_probabilities_options_list : ms_compute_probabilities_option COMMA ms_compute_probabilities_options_list
+                                      | ms_compute_probabilities_option
+                                      ;
+
+ms_compute_probabilities : MS_COMPUTE_PROBABILITIES ';'
+                           { driver.ms_compute_probabilities(); }
+                         | MS_COMPUTE_PROBABILITIES '(' ms_compute_probabilities_options_list ')' ';'
+                           { driver.ms_compute_probabilities(); }
+                         ;
+
+ms_compute_mdd_option : o_output_file_tag
+                      | o_load_mh_file
+                      | o_proposal_type
+                      | o_mdd_proposal_draws
+                      | o_mdd_use_mean_center
                       ;
 
-ms_sbvar : MS_SBVAR ';'
-           { driver.ms_sbvar(); }
-         | MS_SBVAR '(' ms_sbvar_options_list ')' ';'
-           { driver.ms_sbvar(); }
-         ;
+ms_compute_mdd_options_list : ms_compute_mdd_option COMMA ms_compute_mdd_options_list
+                            | ms_compute_mdd_option
+                            ;
+
+ms_compute_mdd : MS_COMPUTE_MDD ';'
+                 { driver.ms_compute_mdd(); }
+               | MS_COMPUTE_MDD '(' ms_compute_mdd_options_list ')' ';'
+                 { driver.ms_compute_mdd(); }
+               ;
+
+ms_simulation_option : o_output_file_tag
+                     | o_mh_replic
+                     | o_drop
+                     | o_thinning_factor
+                     | o_adaptive_mh_draws
+                     ;
 
+ms_simulation_options_list : ms_simulation_option COMMA ms_simulation_options_list
+                           | ms_simulation_option
+                           ;
+
+ms_simulation : MS_SIMULATION ';'
+                { driver.ms_simulation(); }
+              | MS_SIMULATION '(' ms_simulation_options_list ')' ';'
+                { driver.ms_simulation(); }
+              ;
+
+ms_estimation_option : o_coefficients_prior_hyperparameters
+                     | o_freq
+                     | o_initial_year
+                     | o_initial_subperiod
+                     | o_final_year
+                     | o_final_subperiod
+                     | o_datafile
+                     | o_varlist
+                     | o_nlags
+                     | o_cross_restrictions
+                     | o_contemp_reduced_form
+                     | o_bayesian_prior
+                     | o_alpha
+                     | o_beta
+                     | o_gsig2_lmd
+                     | o_gsig2_lmdm
+                     | o_upper_cholesky
+                     | o_lower_cholesky
+                     | o_output_file_tag
+                     | o_convergence_starting_value
+                     | o_convergence_ending_value
+                     | o_convergence_increment_value
+                     | o_max_iterations_starting_value
+                     | o_max_iterations_increment_value
+                     | o_max_block_iterations
+                     | o_max_repeated_optimization_runs
+                     | o_function_convergence_criterion
+                     | o_parameter_convergence_criterion
+                     | o_number_of_large_perturbations
+                     | o_number_of_small_perturbations
+                     | o_number_of_posterior_draws_after_perturbation
+                     | o_max_number_of_stages
+                     | o_random_function_convergence_criterion
+                     | o_random_parameter_convergence_criterion
+                     ;
+
+ms_estimation_options_list : ms_estimation_option COMMA ms_estimation_options_list
+                           | ms_estimation_option
+                           ;
+
+ms_estimation : MS_ESTIMATION ';'
+                { driver.ms_estimation(); }
+              | MS_ESTIMATION '(' ms_estimation_options_list ')' ';'
+                { driver.ms_estimation(); }
+              ;
 
 dynare_sensitivity : DYNARE_SENSITIVITY ';'
                      { driver.dynare_sensitivity(); }
@@ -1896,25 +1998,9 @@ o_eq_cms : EQ_CMS EQUAL INT_NUMBER {driver.option_num("ms.eq_cms",$3); };
 o_tlindx : TLINDX EQUAL INT_NUMBER {driver.option_num("ms.tlindx",$3); };
 o_tlnumber : TLNUMBER EQUAL INT_NUMBER {driver.option_num("ms.tlnumber",$3); };
 o_cnum : CNUM EQUAL INT_NUMBER {driver.option_num("ms.cnum",$3); };
-o_output_file_tag : OUTPUT_FILE_TAG EQUAL '(' symbol_list ')' {driver.option_symbol_list("ms.output_file_tag"); };
-o_create_initialization_file : CREATE_INITIALIZATION_FILE EQUAL INT_NUMBER {driver.option_num("ms.create_initialization_file",$3); };
-o_estimate_msmodel : ESTIMATE_MSMODEL EQUAL INT_NUMBER {driver.option_num("ms.estimate_msmodel",$3); };
-o_compute_mdd : COMPUTE_MDD EQUAL INT_NUMBER {driver.option_num("ms.compute_mdd",$3); };
-o_compute_probabilities : COMPUTE_PROBABILITIES EQUAL INT_NUMBER {driver.option_num("ms.compute_probabilities",$3); };
-o_print_draws : PRINT_DRAWS EQUAL INT_NUMBER {driver.option_num("ms.print_draws",$3); };
-o_n_draws : N_DRAWS EQUAL INT_NUMBER {driver.option_num("ms.n_draws",$3); };
-o_thinning_factor : THINNING_FACTOR EQUAL INT_NUMBER {driver.option_num("ms.thinning_factor",$3); };
-o_markov_file : MARKOV_FILE EQUAL NAME {driver.option_str("ms.markov_file",$3); };
-o_mhm_file: MHM_FILE EQUAL NAME {driver.option_str("ms.mhm_file",$3); };
-o_proposal_draws : PROPOSAL_DRAWS EQUAL INT_NUMBER {driver.option_num("ms.proposal_draws",$3); };
-o_draws_nbr_burn_in_1 : DRAWS_NBR_BURN_IN_1 EQUAL INT_NUMBER {driver.option_num("ms.draws_nbr_burn_in_1",$3); };
-o_draws_nbr_burn_in_2 : DRAWS_NBR_BURN_IN_2 EQUAL INT_NUMBER {driver.option_num("ms.draws_nbr_burn_in_2",$3); };
-o_draws_nbr_mean_var_estimate : DRAWS_NBR_MEAN_VAR_ESTIMATE EQUAL INT_NUMBER {driver.option_num("ms.draws_nbr_mean_var_estimate",$3); };
-o_draws_nbr_modified_harmonic_mean : DRAWS_NBR_MODIFIED_HARMONIC_MEAN EQUAL INT_NUMBER {driver.option_num("ms.draws_nbr_modified_harmonic_mean",$3); };
-o_dirichlet_scale : DIRICHLET_SCALE EQUAL INT_NUMBER {driver.option_num("ms.dirichlet_scale",$3); };
+o_output_file_tag : OUTPUT_FILE_TAG EQUAL '(' symbol_list ')' {driver.option_symbol_list("ms.ft"); };
 o_k_order_solver : K_ORDER_SOLVER {driver.option_num("k_order_solver","1"); };
 o_pruning : PRUNING { driver.option_num("pruning", "1"); };
-
 o_chain : CHAIN EQUAL INT_NUMBER { driver.option_num("ms.chain",$3); };
 o_state : STATE EQUAL INT_NUMBER { driver.option_num("ms.state",$3); };
 o_duration : DURATION EQUAL non_negative_number
@@ -1958,6 +2044,56 @@ o_selected_variables_only : SELECTED_VARIABLES_ONLY
 o_cova_compute : COVA_COMPUTE EQUAL INT_NUMBER
                  { driver.option_num("cova_compute",$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
+                                       { driver.option_num("coefficients_prior_hyperparameters",$3); };
+o_convergence_starting_value : CONVERGENCE_STARTING_VALUE EQUAL non_negative_number
+                               { driver.option_num("ms.convergence_starting_value",$3); };
+o_convergence_ending_value : CONVERGENCE_ENDING_VALUE EQUAL non_negative_number
+                             { driver.option_num("ms.convergence_ending_value",$3); };
+o_convergence_increment_value : CONVERGENCE_INCREMENT_VALUE EQUAL non_negative_number
+                                { driver.option_num("ms.convergence_increment_value",$3); };
+o_max_iterations_starting_value : MAX_ITERATIONS_STARTING_VALUE EQUAL INT_NUMBER
+                                  { driver.option_num("ms.max_iterations_starting_value",$3); };
+o_max_iterations_increment_value : MAX_ITERATIONS_INCREMENT_VALUE EQUAL non_negative_number
+                                   { driver.option_num("ms.max_iterations_increment_value",$3); };
+o_max_block_iterations : MAX_BLOCK_ITERATIONS EQUAL INT_NUMBER
+                         { driver.option_num("ms.max_block_iterations",$3); };
+o_max_repeated_optimization_runs : MAX_REPEATED_OPTIMIZATION_RUNS EQUAL INT_NUMBER
+                                   { driver.option_num("ms.max_repeated_optimization_runs",$3); };
+o_function_convergence_criterion : FUNCTION_CONVERGENCE_CRITERION EQUAL non_negative_number
+                                   { driver.option_num("ms.function_convergence_criterion",$3); };
+o_parameter_convergence_criterion : PARAMETER_CONVERGENCE_CRITERION EQUAL non_negative_number
+                                    { driver.option_num("ms.parameter_convergence_criterion",$3); };
+o_number_of_large_perturbations : NUMBER_OF_LARGE_PERTURBATIONS EQUAL INT_NUMBER
+                                  { driver.option_num("ms.number_of_large_perturbations",$3); };
+o_number_of_small_perturbations : NUMBER_OF_SMALL_PERTURBATIONS EQUAL INT_NUMBER
+                                  { driver.option_num("ms.number_of_small_perturbations",$3); };
+o_number_of_posterior_draws_after_perturbation : NUMBER_OF_POSTERIOR_DRAWS_AFTER_PERTURBATION EQUAL INT_NUMBER
+                                                 { driver.option_num("ms.number_of_posterior_draws_after_perturbation",$3); };
+o_max_number_of_stages : MAX_NUMBER_OF_STAGES EQUAL INT_NUMBER
+                         { driver.option_num("ms.max_number_of_stages",$3); };
+o_random_function_convergence_criterion : RANDOM_FUNCTION_CONVERGENCE_CRITERION EQUAL non_negative_number
+                                          { driver.option_num("ms.random_function_convergence_criterion",$3); };
+o_random_parameter_convergence_criterion : RANDOM_PARAMETER_CONVERGENCE_CRITERION EQUAL non_negative_number
+                                           { driver.option_num("ms.random_parameter_convergence_criterion",$3); };
+o_thinning_factor : THINNING_FACTOR EQUAL INT_NUMBER { driver.option_num("ms.thinning_factor",$3); };
+o_adaptive_mh_draws : ADAPTIVE_MH_DRAWS EQUAL INT_NUMBER { driver.option_num("ms.adaptive_mh_draws",$3); };
+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_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_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"); };
 
 range : symbol ':' symbol
         {
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index e27b8ee92a6be042a3bc2348ca72165fcc8ababd..d1a4b452de921ecd198275fa84258294d4bbba8d 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -140,7 +140,13 @@ string eofbuff;
 <INITIAL>forecast {BEGIN DYNARE_STATEMENT; return token::FORECAST;}
 <INITIAL>shock_decomposition {BEGIN DYNARE_STATEMENT; return token::SHOCK_DECOMPOSITION;}
 <INITIAL>sbvar {BEGIN DYNARE_STATEMENT; return token::SBVAR;}
-<INITIAL>ms_sbvar {BEGIN DYNARE_STATEMENT; return token::MS_SBVAR;}
+<INITIAL>ms_estimation {BEGIN DYNARE_STATEMENT; return token::MS_ESTIMATION;}
+<INITIAL>ms_simulation {BEGIN DYNARE_STATEMENT; return token::MS_SIMULATION;}
+<INITIAL>ms_compute_mdd {BEGIN DYNARE_STATEMENT; return token::MS_COMPUTE_MDD;}
+<INITIAL>ms_compute_probabilities {BEGIN DYNARE_STATEMENT; return token::MS_COMPUTE_PROBABILITIES;}
+<INITIAL>ms_forecast {BEGIN DYNARE_STATEMENT; return token::MS_FORECAST;}
+<INITIAL>ms_irf {BEGIN DYNARE_STATEMENT; return token::MS_IRF;}
+<INITIAL>ms_variance_decomposition {BEGIN DYNARE_STATEMENT; return token::MS_VARIANCE_DECOMPOSITION;}
 <INITIAL>conditional_forecast {BEGIN DYNARE_STATEMENT; return token::CONDITIONAL_FORECAST;}
 <INITIAL>plot_conditional_forecast {BEGIN DYNARE_STATEMENT; return token::PLOT_CONDITIONAL_FORECAST;}
 
@@ -312,22 +318,41 @@ string eofbuff;
   return token::CNUM;
 }
 <DYNARE_STATEMENT>banact {return token::BANACT;}
+
 <DYNARE_STATEMENT>output_file_tag {return token::OUTPUT_FILE_TAG;}
-<DYNARE_STATEMENT>create_initialization_file {return token::CREATE_INITIALIZATION_FILE;}
-<DYNARE_STATEMENT>estimate_msmodel {return token::ESTIMATE_MSMODEL;}
-<DYNARE_STATEMENT>compute_mdd {return token::COMPUTE_MDD;}
-<DYNARE_STATEMENT>compute_probabilities {return token::COMPUTE_PROBABILITIES;}
-<DYNARE_STATEMENT>print_draws {return token::PRINT_DRAWS;}
-<DYNARE_STATEMENT>n_draws {return token::N_DRAWS;}
+<DYNARE_STATEMENT>filtered {return token::FILTERED;}
+<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>thinning_factor {return token::THINNING_FACTOR;}
-<DYNARE_STATEMENT>markov_file {return token::MARKOV_FILE;}
-<DYNARE_STATEMENT>mhm_file {return token::MHM_FILE;}
-<DYNARE_STATEMENT>proposal_draws {return token::PROPOSAL_DRAWS;}
-<DYNARE_STATEMENT>draws_nbr_burn_in_1 {return token::DRAWS_NBR_BURN_IN_1;}
-<DYNARE_STATEMENT>draws_nbr_burn_in_2 {return token::DRAWS_NBR_BURN_IN_2;}
-<DYNARE_STATEMENT>draws_nbr_mean_var_estimate {return token::DRAWS_NBR_MEAN_VAR_ESTIMATE;}
-<DYNARE_STATEMENT>draws_nbr_modified_harmonic_mean {return token::DRAWS_NBR_MODIFIED_HARMONIC_MEAN;}
-<DYNARE_STATEMENT>dirichlet_scale {return token::DIRICHLET_SCALE;}
+<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;}
+<DYNARE_STATEMENT>filtered_probabilities {return token::FILTERED_PROBABILITIES;}
+<DYNARE_STATEMENT>real_time_smoothed {return token::REAL_TIME_SMOOTHED;}
+<DYNARE_STATEMENT>proposal_type {return token::PROPOSAL_TYPE;}
+<DYNARE_STATEMENT>mdd_proposal_draws {return token::MDD_PROPOSAL_DRAWS;}
+<DYNARE_STATEMENT>mdd_use_mean_center {return token::MDD_USE_MEAN_CENTER;}
+<DYNARE_STATEMENT>adaptive_mh_draws {return token::ADAPTIVE_MH_DRAWS;}
+<DYNARE_STATEMENT>coefficients_prior_hyperparameters {return token::COEFFICIENTS_PRIOR_HYPERPARAMETERS;}
+<DYNARE_STATEMENT>convergence_starting_value {return token::CONVERGENCE_STARTING_VALUE;}
+<DYNARE_STATEMENT>convergence_ending_value {return token::CONVERGENCE_ENDING_VALUE;}
+<DYNARE_STATEMENT>convergence_increment_value {return token::CONVERGENCE_INCREMENT_VALUE;}
+<DYNARE_STATEMENT>max_iterations_starting_value {return token::MAX_ITERATIONS_STARTING_VALUE;}
+<DYNARE_STATEMENT>max_iterations_increment_value {return token::MAX_ITERATIONS_INCREMENT_VALUE;}
+<DYNARE_STATEMENT>max_block_iterations {return token::MAX_BLOCK_ITERATIONS;}
+<DYNARE_STATEMENT>max_repeated_optimization_runs {return token::MAX_REPEATED_OPTIMIZATION_RUNS;}
+<DYNARE_STATEMENT>function_convergence_criterion {return token::FUNCTION_CONVERGENCE_CRITERION;}
+<DYNARE_STATEMENT>parameter_convergence_criterion {return token::PARAMETER_CONVERGENCE_CRITERION;}
+<DYNARE_STATEMENT>number_of_large_perturbations {return token::NUMBER_OF_LARGE_PERTURBATIONS;}
+<DYNARE_STATEMENT>number_of_small_perturbations {return token::NUMBER_OF_SMALL_PERTURBATIONS;}
+<DYNARE_STATEMENT>number_of_posterior_draws_after_perturbation {return token::NUMBER_OF_POSTERIOR_DRAWS_AFTER_PERTURBATION;}
+<DYNARE_STATEMENT>max_number_of_stages {return token::MAX_NUMBER_OF_STAGES;}
+<DYNARE_STATEMENT>random_function_convergence_criterion {return token::RANDOM_FUNCTION_CONVERGENCE_CRITERION;}
+<DYNARE_STATEMENT>random_parameter_convergence_criterion {return token::RANDOM_PARAMETER_CONVERGENCE_CRITERION;}
+
 <DYNARE_STATEMENT>instruments {return token::INSTRUMENTS;}
 
  /* These four (var, varexo, varexo_det, parameters) are for change_type */
@@ -419,8 +444,8 @@ string eofbuff;
 <DYNARE_BLOCK>equation {return token::EQUATION;}
 <DYNARE_BLOCK>exclusion {return token::EXCLUSION;}
 <DYNARE_BLOCK>lag {return token::LAG;}
-<DYNARE_BLOCK>upper_cholesky {return token::UPPER_CHOLESKY;}
-<DYNARE_BLOCK>lower_cholesky {return token::LOWER_CHOLESKY;}
+<DYNARE_STATEMENT,DYNARE_BLOCK>upper_cholesky {return token::UPPER_CHOLESKY;}
+<DYNARE_STATEMENT,DYNARE_BLOCK>lower_cholesky {return token::LOWER_CHOLESKY;}
 <DYNARE_STATEMENT>chain {return token::CHAIN;}
 <DYNARE_STATEMENT>state {return token::STATE;}
 <DYNARE_STATEMENT>number_of_states {return token::NUMBER_OF_STATES;}
diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc
index ac99d35683272b5eefa2d36a27704adec2fc9e4b..e624b8198325e99af63a15ff1dc32af1be8f7404 100644
--- a/preprocessor/ParsingDriver.cc
+++ b/preprocessor/ParsingDriver.cc
@@ -1312,9 +1312,51 @@ ParsingDriver::sbvar()
 }
 
 void
-ParsingDriver::ms_sbvar()
+ParsingDriver::ms_estimation()
 {
-  mod_file->addStatement(new MS_SBVARStatement(options_list));
+  mod_file->addStatement(new MSSBVAREstimationStatement(options_list));
+  options_list.clear();
+}
+
+void
+ParsingDriver::ms_simulation()
+{
+  mod_file->addStatement(new MSSBVARSimulationStatement(options_list));
+  options_list.clear();
+}
+
+void
+ParsingDriver::ms_compute_mdd()
+{
+  mod_file->addStatement(new MSSBVARComputeMDDStatement(options_list));
+  options_list.clear();
+}
+
+void
+ParsingDriver::ms_compute_probabilities()
+{
+  mod_file->addStatement(new MSSBVARComputeProbabilitiesStatement(options_list));
+  options_list.clear();
+}
+
+void
+ParsingDriver::ms_irf()
+{
+  mod_file->addStatement(new MSSBVARIrfStatement(options_list));
+  options_list.clear();
+}
+
+void
+ParsingDriver::ms_forecast()
+{
+  mod_file->addStatement(new MSSBVARForecastStatement(options_list));
+  options_list.clear();
+}
+
+void
+ParsingDriver::ms_variance_decomposition()
+{
+  mod_file->addStatement(new MSSBVARVarianceDecompositionStatement(options_list));
   options_list.clear();
 }
 
diff --git a/preprocessor/ParsingDriver.hh b/preprocessor/ParsingDriver.hh
index 3140ada7f6872afb2f1f5a2f8c5e6820f155ff55..990733917dce2e141595201505d299bf2c3880d9 100644
--- a/preprocessor/ParsingDriver.hh
+++ b/preprocessor/ParsingDriver.hh
@@ -392,8 +392,20 @@ public:
   void bvar_forecast(string *nlags);
   //! SBVAR statement
   void sbvar();
-  //! MS_SBVAR statement
-  void ms_sbvar();
+  //! Markov Switching Statement: Estimation
+  void ms_estimation();
+  //! Markov Switching Statement: Simulation
+  void ms_simulation();
+  //! Markov Switching Statement: MDD
+  void ms_compute_mdd();
+  //! Markov Switching Statement: Probabilities
+  void ms_compute_probabilities();
+  //! Markov Switching Statement: IRF
+  void ms_irf();
+  //! Markov Switching Statement: Forecast
+  void ms_forecast();
+  //! Markov Switching Statement: Variance Decomposition
+  void ms_variance_decomposition();
   //! Svar statement
   void svar();
   //! MarkovSwitching statement
diff --git a/tests/ms-sbvar/test_ms_variances.mod b/tests/ms-sbvar/test_ms_variances.mod
index 18303862451109b59a2f98dae936c5043edff2d2..1380683def1f7adc1410ac44ef810686b148007a 100644
--- a/tests/ms-sbvar/test_ms_variances.mod
+++ b/tests/ms-sbvar/test_ms_variances.mod
@@ -1,4 +1,3 @@
-// same as test_upper_cholesky.mod, but with reordered variables. Results must be the same.
 var R Pie Y;
 
 model;
@@ -17,5 +16,7 @@ markov_switching(chain=1,number_of_states=2,duration=2.5);
 
 svar(variances, chain=1);
 
-ms_sbvar(datafile = data,freq=4,initial_year=1959,final_year=2005,nlags=4,draws_nbr_modified_harmonic_mean=10000);
-
+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;