diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 8b8029c7f..207115da8 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -283,6 +283,15 @@ options_.ms.print_draws = 1;
 options_.ms.n_draws=1000;
 options_.ms.thinning_factor=1;
 options_.ms.proposal_draws = 100000;
-
+options_.ms.lower_cholesky = 0;
+options_.ms.upper_cholesky = 0;
+options_.ms.Qi = [];
+options_.ms.Ri = [];
+options_.ms.draws_nbr_burn_in_1 = 30000;
+options_.ms.draws_nbr_burn_in_2 = 10000;
+options_.ms.draws_nbr_mean_var_estimate = 200000;
+options_.ms.draws_nbr_modified_harmonic_mean = 1000000;
+options_.ms.thinning_factor = 1;
+options_.ms.dirichlet_scale = 1;
 % initialize persistent variables in priordens()
 priordens([],[],[],[],[],[],1);
diff --git a/matlab/swz/bin/.gitgnore b/matlab/swz/bin/.gitgnore
new file mode 100644
index 000000000..e69de29bb
diff --git a/matlab/swz/c-code/Makefile b/matlab/swz/c-code/Makefile
index a0fb833f4..b39b51f5b 100644
--- a/matlab/swz/c-code/Makefile
+++ b/matlab/swz/c-code/Makefile
@@ -50,10 +50,11 @@ endif
 
 ifdef USE_MICHEL_LAPTOP
 
-#CC = gcc
-CC = /opt/intel/Compiler/11.0/074/bin/intel64/icc -no-multibyte-chars
+CC = gcc
+#CC = /opt/intel/Compiler/11.0/074/bin/intel64/icc -no-multibyte-chars
 FC = gfortran
 CFLAGS = -g -static
+#CFLAGS = -g
 
 TAO_DIR = $(WORK_DIR)/utilities/TZCcode
 
@@ -65,6 +66,7 @@ INTEL_LIBS = -lguide
 MKL_BASE_DIR = /opt/intel/Compiler/11.0/074/mkl
 MKL_LIBS_DIR = $(MKL_BASE_DIR)/lib/em64t
 MKL_LIBS = -lmkl_lapack -lmkl_em64t 
+
 endif
 
 
@@ -173,15 +175,16 @@ ERROR_DIR = $(WORK_DIR)/utilities/DWCcode/error
 ARRAY_DIR = $(WORK_DIR)/utilities/DWCcode/arrays
 ASCII_DIR = $(WORK_DIR)/utilities/DWCcode/ascii 
 STAT_DIR = $(WORK_DIR)/utilities/DWCcode/stat
+SPHERICAL_DIR = $(WORK_DIR)/utilities/DWCcode/spherical
 SORT_DIR = $(WORK_DIR)/utilities/DWCcode/sort
 SWITCH_DIR = $(WORK_DIR)/sbvar/switching
 VAR_DIR = $(WORK_DIR)/sbvar/var
 #################################################################################
 
 # DW  FILES
-INCLUDE_DIR := $(INCLUDE_DIR) -I$(MATRIX_DIR) -I$(ERROR_DIR) -I$(ARRAY_DIR) -I$(ASCII_DIR) -I$(STAT_DIR) -I$(SWITCH_DIR) -I$(VAR_DIR)
-VPATH := $(VPATH) $(MATRIX_DIR) $(ERROR_DIR) $(ARRAY_DIR) $(ASCII_DIR) $(STAT_DIR) $(SWITCH_DIR) $(VAR_DIR)
-OBJS := $(OBJS) bmatrix.o matrix.o dw_error.o dw_rand.o dw_matrix_rand.o dw_array.o dw_matrix_array.o dw_ascii.o dw_parse_cmd.o 
+INCLUDE_DIR := -I$(MATRIX_DIR) -I$(ERROR_DIR) -I$(ARRAY_DIR) -I$(ASCII_DIR) -I$(STAT_DIR) -I$(SPHERICAL_DIR) -I$(SORT_DIR) -I$(SWITCH_DIR) -I$(VAR_DIR) $(INCLUDE_DIR) 
+VPATH := $(VPATH) $(MATRIX_DIR) $(ERROR_DIR) $(ARRAY_DIR) $(ASCII_DIR) $(STAT_DIR) $(SPHERICAL_DIR) $(SORT_DIR) $(SWITCH_DIR) $(VAR_DIR)
+OBJS := $(OBJS) bmatrix.o matrix.o dw_error.o dw_rand.o dw_matrix_rand.o dw_array.o dw_matrix_array.o dw_matrix_sort.o dw_ascii.o dw_parse_cmd.o
 
 # TAO FILES
 OBJS := $(OBJS)
@@ -189,29 +192,38 @@ OBJS := $(OBJS)
 # PROJECT FILE
 INCLUDE_DIR := $(INCLUDE_DIR)
 VPATH := $(VPATH)
-OBJS := $(OBJS) PrintDraws.o switch.o switchio.o VARbase.o VARio.o command_line_VAR.o
+OBJS_DRAWS := $(OBJS) PrintDraws.o switch.o switchio.o VARbase.o VARio.o command_line_VAR.o
+OBJS_TAO := tzmatlab.o mathlib.o cstz_dw.o
+OBJS_ESTIM := $(OBJS) $(OBJS_TAO) estimate.o VARbase.o VARio.o command_line_VAR.o switch.o switchio.o switch_opt.o csminwel.o
+OBJS_INIT := $(OBJS) create_init_file.o switch.o switchio.o VARbase.o VARio.o VARio_matlab.o 
+OBJS_MHM_1 := $(OBJS)  mhm_VAR_main_1.o mhm_VAR.o VARbase.o VARio.o command_line_VAR.o switch.o switchio.o 
+OBJS_MHM_2 := $(OBJS) mhm_VAR_main_2.o spherical.o VARbase.o VARio.o switch.o switchio.o mhm_VAR.o
+OBJS_PROBA := $(OBJS) probabilities.o switch.o switchio.o VARbase.o VARio.o command_line_VAR.o
+
+
 
 #################################################################################
 # OUTPUT
 
-all: $(OUT_DIR)/sbvar_draws $(OUT_DIR)/sbvar_estimation $(OUT_DIR)/sbvar_init_file $(OUT_DIR)/sbvar_mhm_1 $(OUT_DIR)/sbvar_mhm_2 $(OUT_DIR)/sbvar_probabilities
+#all: $(OUT_DIR)/sbvar_draws $(OUT_DIR)/sbvar_estimation $(OUT_DIR)/sbvar_init_file $(OUT_DIR)/sbvar_mhm_1 $(OUT_DIR)/sbvar_mhm_2 $(OUT_DIR)/sbvar_probabilities
+all: $(OUT_DIR)/sbvar_draws $(OUT_DIR)/sbvar_init_file $(OUT_DIR)/sbvar_mhm_1 $(OUT_DIR)/sbvar_mhm_2 $(OUT_DIR)/sbvar_probabilities
 
-$(OUT_DIR)/sbvar_draws: $(OBJS)
+$(OUT_DIR)/sbvar_draws: $(OBJS_DRAWS)
 	$(CC) $(CFLAGS) $^ $(LIBS_DIR) $(LIBS) -o $(OUT_DIR)/sbvar_draws
 
-$(OUT_DIR)/sbvar_estimation: $(OBJS)
+$(OUT_DIR)/sbvar_estimation: $(OBJS_ESTIM)
 	$(CC) $(CFLAGS) $^ $(LIBS_DIR) $(LIBS) -o $(OUT_DIR)/sbvar_estimation
 
-$(OUT_DIR)/sbvar_init_file: $(OBJS)
+$(OUT_DIR)/sbvar_init_file: $(OBJS_INIT)
 	$(CC) $(CFLAGS) $^ $(LIBS_DIR) $(LIBS) -o $(OUT_DIR)/sbvar_init_file
 
-$(OUT_DIR)/sbvar_mhm_1: $(OBJS)
+$(OUT_DIR)/sbvar_mhm_1: $(OBJS_MHM_1)
 	$(CC) $(CFLAGS) $^ $(LIBS_DIR) $(LIBS) -o $(OUT_DIR)/sbvar_mhm_1
 
-$(OUT_DIR)/sbvar_mhm_2: $(OBJS)
+$(OUT_DIR)/sbvar_mhm_2: $(OBJS_MHM_2)
 	$(CC) $(CFLAGS) $^ $(LIBS_DIR) $(LIBS) -o $(OUT_DIR)/sbvar_mhm_2
 
-$(OUT_DIR)/sbvar_probabilities: $(OBJS)
+$(OUT_DIR)/sbvar_probabilities: $(OBJS_PROBA)
 	$(CC) $(CFLAGS) $^ $(LIBS_DIR) $(LIBS) -o $(OUT_DIR)/sbvar_probabilities
 
 
diff --git a/matlab/swz/identification/lower_cholesky.m b/matlab/swz/identification/lower_cholesky.m
new file mode 100644
index 000000000..0ce27ce5f
--- /dev/null
+++ b/matlab/swz/identification/lower_cholesky.m
@@ -0,0 +1,32 @@
+function [Ui,Vi,n0,np,ixmC0Pres] = upper_cholesky(nvar,nexo,options_ms)
+
+lags = options_ms.nlags;
+indxC0Pres = options_ms.cross_restrictions;
+
+Ui = cell(nvar,1);
+Vi = cell(nvar,1);
+n0 = zeros(nvar,1);
+np = zeros(nvar,1);
+if (nargin==2)
+   nexo = 1; 
+elseif (nargin==3)
+   indxC0Pres = 0;
+end
+k = lags*nvar+nexo; 
+Qi = zeros(nvar,nvar,nvar);
+Ri = zeros(k,k,nvar);
+
+for ii=2:nvar
+    Pi=diag(diag(ones(nvar-ii+1)),ii-1);
+    Qi(:,:,ii-1)=Pi(1:nvar,1:nvar);
+    Qi(:,:,nvar)=zeros(nvar,nvar);
+end
+
+for n=1:nvar 
+   Ui{n} = null(Qi(:,:,n));
+   Vi{n} = null(Ri(:,:,n));
+   n0(n) = size(Ui{n},2);
+   np(n) = size(Vi{n},2);
+end
+
+ixmC0Pres = NaN;
\ No newline at end of file
diff --git a/matlab/swz/identification/upper_cholesky.m b/matlab/swz/identification/upper_cholesky.m
index 8c41339ac..471d80a7f 100644
--- a/matlab/swz/identification/upper_cholesky.m
+++ b/matlab/swz/identification/upper_cholesky.m
@@ -1,22 +1,28 @@
-function [Ui,Vi,n0,np,ixmC0Pres] = upper_cholesky(lags,nvar,nexo,indxC0Pres)
+function [Ui,Vi,n0,np,ixmC0Pres] = lower_cholesky(nvar,nexo,options_ms)
+
+lags = options_ms.nlags;
+indxC0Pres = options_ms.cross_restrictions;
 
 Ui = cell(nvar,1);
 Vi = cell(nvar,1);
 n0 = zeros(nvar,1);
 np = zeros(nvar,1);
+
 if (nargin==2)
-   nexo = 1; 
+   nexo = 1;
 elseif (nargin==3)
    indxC0Pres = 0;
 end
-k = lags*nvar+nexo; 
+
+k = lags*nvar+nexo;
 Qi = zeros(nvar,nvar,nvar);
 Ri = zeros(k,k,nvar);
 
 for ii=2:nvar
-    Pi=diag(diag(ones(nvar-ii+1)),ii-1);
-    Qi(:,:,ii-1)=Pi(1:nvar,1:nvar);
-    Qi(:,:,nvar)=zeros(nvar,nvar);
+    
+  Qi(ii-1,ii-1,ii)=1;
+  Qi(:,:,ii)=Qi(:,:,ii)+Qi(:,:,ii-1);
+   
 end
 
 for n=1:nvar 
diff --git a/matlab/swz/msstart2.m b/matlab/swz/msstart2.m
index 31012b85b..b72be117e 100644
--- a/matlab/swz/msstart2.m
+++ b/matlab/swz/msstart2.m
@@ -71,7 +71,7 @@ IndxNmlr = [1 0 0 0 0 0];  % imported by nmlzvar.m
 %
 %options_.ms.indxgdls = 1;  % 1: graph point forecast on variables; 0: disable
 nconstr1=nfqm;      % number of the 1st set of constraints
-nconstr2=options_.ms.forecast ;     % number of the 2nd set of constraints
+nconstr2=options_.forecast ;     % number of the 2nd set of constraints
 nconstr=nconstr1+nconstr2;   % q: 4 years -- 4*12 months.
                          % When 0, no conditions directly on variables <<>>
 nconstr=0   %6*nconstr1;
@@ -274,7 +274,7 @@ if options_.ms.indxestima
    ye = [dateswd y];
 
    %* Obtain linear restrictions
-   [Uiconst,Viconst,n0,np,ixmC0Pres] = feval(options_.ms.restriction_fname,options_.ms.nlags ,nvar,nexo,options_.ms.cross_restrictions );
+   [Uiconst,Viconst,n0,np,ixmC0Pres] = feval(options_.ms.restriction_fname,nvar,nexo,options_.ms );
    if min(n0)==0
       disp(' ')
       warning('A0: restrictions in dlrprior.m give no free parameter in one of equations')
@@ -574,7 +574,7 @@ if options_.ms.indxestima
 
          %*** conditions in every period
          vpntM = fn_dataext(E1yrqm, E2yrqm,xdatae); % point value matrix with dates
-         %     vaveM = fn_dataext([yrEnd+1 0],[yrEnd+options_.ms.forecast  0],yact2yre); % average value matrix with dates
+         %     vaveM = fn_dataext([yrEnd+1 0],[yrEnd+options_.forecast  0],yact2yre); % average value matrix with dates
          for i=1:nconstr
             if i<=nconstr1
                valuecon(i) = vpntM(i,2+varcon(i)); % 2: first 2 elements are dates
@@ -593,10 +593,10 @@ if options_.ms.indxestima
 
 %         %*** average condtions over,say, options_.ms.freq periods.
 %         if qmEnd==options_.ms.freq
-%            vaveM = fn_dataext([yrEnd+1 0],[yrEnd+options_.ms.forecast  0],yact2yre); % average value matrix with dates
+%            vaveM = fn_dataext([yrEnd+1 0],[yrEnd+options_.forecast  0],yact2yre); % average value matrix with dates
 %            valuecon(1) = vaveM(1,2+varcon(1));  % 2: first 2 elements are dates
 %         else
-%            vaveM = fn_dataext([yrEnd 0],[yrEnd+options_.ms.forecast  0],yact2yre); % average value matrix with dates
+%            vaveM = fn_dataext([yrEnd 0],[yrEnd+options_.forecast  0],yact2yre); % average value matrix with dates
 %            yactrem = fn_dataext([yrEnd qmEnd+1],[yrEnd options_.ms.freq],xdatae);
 %            valuecon(1) = sum(yactrem(:,2+varcon(1)),1)/length(stepcon{1});
 %                                    % 2: first 2 elements are dates
diff --git a/matlab/swz/msstart_setup.m b/matlab/swz/msstart_setup.m
index 6411e0701..e41559602 100644
--- a/matlab/swz/msstart_setup.m
+++ b/matlab/swz/msstart_setup.m
@@ -69,8 +69,8 @@ yrStart=options_.ms.initial_year;
 qmStart=options_.ms.initial_subperiod;
 yrEnd=options_.ms.final_year;
 qmEnd=options_.ms.final_subperiod;
-%options_.ms.forecast = 4;   % number of years for forecasting
-if options_.ms.forecast<1
+%options_.forecast = 4;   % number of years for forecasting
+if options_.forecast<1
    error('To be safe, the number of forecast years should be at least 1')
 end
 ystr=num2str(yrEnd);
@@ -82,7 +82,7 @@ if qmEnd==options_.ms.freq
 else
    E1yrqm = [yrEnd qmEnd+1];  % first year and quarter (month) after the sample
 end
-E2yrqm = [yrEnd+options_.ms.forecast qmEnd];   % end at the last month (quarter) of a calendar year after the sample
+E2yrqm = [yrEnd+options_.forecast qmEnd];   % end at the last month (quarter) of a calendar year after the sample
 [fdates,nfqm]=fn_calyrqm(options_.ms.freq,E1yrqm,E2yrqm);   % forecast dates and number of forecast dates
 [sdates,nsqm] = fn_calyrqm(options_.ms.freq,[yrStart qmStart],[yrEnd qmEnd]);
    % sdates: dates for the whole sample (including options_.ms.nlags)
@@ -92,7 +92,7 @@ if nSample~=nsqm
    pause
 end
 imstp = 4*options_.ms.freq;    % <<>>  impulse responses (4 years)
-nayr = 4; %options_.ms.forecast;  % number of years before forecasting for plotting.
+nayr = 4; %options_.forecast;  % number of years before forecasting for plotting.
 
 
 %------- Prior, etc. -------
@@ -132,7 +132,7 @@ indxmsmdeqn = [0; 0; 0; 0];  %This option disenable using this in fn_rnrprior_co
 
 
 tdf = 3;          % degrees of freedom for t-dist for initial draw of the MC loop
-nbuffer = 100;        % a block or buffer of draws (buffer) that is saved to the disk (not memory)
+nbuffer = 1000;        % a block or buffer of draws (buffer) that is saved to the disk (not memory)
 ndraws1=1*nbuffer;         % 1st part of Monte Carlo draws
 ndraws2=10*ndraws1         % 2nd part of Monte Carlo draws
 seednumber = 0; %7910;    %472534;   % if 0, random state at each clock time
diff --git a/matlab/swz/swz_sbvar.m b/matlab/swz/swz_sbvar.m
index bdea69b89..c501fb462 100644
--- a/matlab/swz/swz_sbvar.m
+++ b/matlab/swz/swz_sbvar.m
@@ -10,14 +10,43 @@ addpath([swz_root '/mhm_specification']);
 
 options.data = read_variables(options.datafile,options.varobs,[],options.xls_sheet,options.xls_range);
 
+if options.forecast == 0
+    options.forecast = 4;
+end
+
 options.ms.output_file_tag = M.fname;
 %options.ms.markov_file = 'specification_2v2c.dat';
 %options.ms.mhm_file = 'MHM_input.dat';
 %options.ms.restriction_fname = 'ftd_upperchol3v'; 
 
+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';
+end
 
 if ms_flag == 1
-    sz_prd(options)
+    % changing some option names to match SWZ code 
+    options.ms.firstMetrop = options.ms.draws_nbr_burn_in_1;
+    options.ms.secondMetrop = options.ms.draws_nbr_burn_in_2;
+    options.ms.ndrawsmv = options.ms.draws_nbr_mean_var_estimate;
+    options.ms.ndrawsmhm = options.ms.draws_nbr_modified_harmonic_mean;
+    options.ms.tfmhm = options.ms.thinning_factor;
+    options.ms.svd = options.ms.dirichlet_scale;
+    % are these options necessary ?
+    options.ms.opt1 = 1;
+    options.ms.opt2 = 1;
+    options.ms.opt3 = 1;
+    % temporary fix
+    options.ms.markov_file = 'markov_file';
+    sz_prd(M,options);
 else
     swz_mardd(options);
 end
diff --git a/matlab/swz/swz_write_markov_file.m b/matlab/swz/swz_write_markov_file.m
new file mode 100644
index 000000000..2a5b64f09
--- /dev/null
+++ b/matlab/swz/swz_write_markov_file.m
@@ -0,0 +1,126 @@
+function swz_write_markov_file(fname,M,options)
+
+    n_chains = length(options.ms.ms_chain);
+    nvars = size(options.varobs,1);
+    
+    fh = fopen(fname,'w');
+    %/******************************************************************************/
+    %/********************* Markov State Variable Information **********************/
+    %/******************************************************************************/
+    
+    fprintf(fh,'//== Flat Independent Markov States and Simple Restrictions ==//\n\n');
+    
+    
+    %//This number is NOT used but read in.
+    fprintf(fh,'//== Number Observations ==//\n');
+    fprintf(fh,'0\n\n');
+
+    fprintf(fh,'//== Number Independent State Variables ==//\n');
+    fprintf(fh,'%d\n\n',n_chains);
+
+    for i_chain = 1:n_chains
+
+        %//=====================================================//
+        %//== state_variable[i] (1 <= i <= n_state_variables) ==//
+        %//=====================================================//
+        fprintf(fh,'//== Number of states for state_variable[%d] ==//\n', ...
+                i_chain);
+        n_states = length(options.ms.ms_chain(i_chain).state);
+        fprintf(fh,'%d\n\n',n_states);
+
+        %//== 03/15/06: DW TVBVAR code reads the data below and overwrite the prior data read somewhere else if any.
+        %//== Each column contains the parameters for a Dirichlet prior on the corresponding
+        %//== column of the transition matrix.  Each element must be positive.  For each column,
+        %//== the relative size of the prior elements determine the relative size of the elements
+        %//== of the transition matrix and overall larger sizes implies a tighter prior.
+        fprintf(fh,['//== Transition matrix prior for state_variable[%d]. ' ...
+                    '(n_states x n_states) ==//\n'],i_chain);
+        Alpha = ones(n_states,n_states);
+        for i_state = 1:n_states
+            p = 1-1/options.ms.ms_chain(i_chain).state(i_state).duration;
+            Alpha(i_state,i_state) = p*(n_states-1)/(1-p);
+            fprintf(fh,'%22.16f',Alpha(i_state,:));
+            fprintf(fh,'\n');
+        end
+ 
+        fprintf(fh,['\n//== Free Dirichet dimensions for state_variable[%d]  ' ...
+                    '==//\n'],i_chain);
+        fprintf(fh,'%d ',repmat(n_states,1,n_states));
+        fprintf(fh,'\n\n');
+
+        %//== The jth restriction matrix is n_states-by-free[j].  Each row of the restriction
+        %//== matrix has exactly one non-zero entry and the sum of each column must be one.
+        fprintf(fh,['//== Column restrictions for state_variable[%d] ' ...
+                    '==//\n'],i_chain);
+        if n_states == 2
+            M = eye(2);
+            for i_state = 1:2
+                for j_state = 1:2
+                    fprintf(fh,'%d ',M(j_state,:));
+                    fprintf(fh,'\n');
+                end
+                fprintf(fh,'\n');
+            end
+        else
+        end
+    end
+
+    %/******************************************************************************/
+    %/******************************* VAR Parameters *******************************/
+    %/******************************************************************************/
+    %//NOT read
+    fprintf(fh,'//== Number Variables ==//\n');
+    fprintf(fh,'%d\n\n',nvars);
+
+    %//NOT read
+    fprintf(fh,'//== Number Lags ==//\n');
+    fprintf(fh,'%d\n\n',options.ms.nlags);
+
+    %//NOT read
+    fprintf(fh,'//== Exogenous Variables ==//\n');
+    fprintf(fh,'1\n\n');
+
+
+    %//== nvar x n_state_variables matrix.  In the jth row, a non-zero value implies that
+    %this state variable controls the jth column of A0 and Aplus
+    fprintf(fh,['//== Controlling states variables for coefficients ==//\' ...
+                'n']);
+    
+    for i_var = 1:nvars
+        for i_chain = 1:n_chains
+            if ~isfield(options.ms.ms_chain(i_chain),'svar_coefficients')
+                i_equations = 0;
+            else
+                i_equations = ...
+                    options.ms.ms_chain(i_chain).svar_coefficients.equations; 
+            end
+            if strcmp(i_equations,'ALL') || any(i_equations == i_var)
+                fprintf(fh,'%d ',1);
+            else
+                fprintf(fh,'%d ',0);
+            end
+        end
+        fprintf(fh,'\n');
+    end
+
+    %//== nvar x n_state_variables matrix.  In the jth row, a non-zero value implies that
+    %this state variable controls the jth diagonal element of Xi
+    fprintf(fh,'\n//== Controlling states variables for variance ==//\n');
+    for i_var = 1:nvars
+        for i_chain = 1:n_chains
+            if ~isfield(options.ms.ms_chain(i_chain),'svar_variances')
+                i_equations = 0;
+            else
+                i_equations = ...
+                    options.ms.ms_chain(i_chain).svar_variances.equations; 
+            end
+            if strcmp(i_equations,'ALL') || any(i_equations == i_var)
+                fprintf(fh,'%d ',1);
+            else
+                fprintf(fh,'%d ',0);
+            end
+        end
+        fprintf(fh,'\n');
+    end
+
+    fclose(fh);
diff --git a/matlab/swz/swz_write_mhm_input.m b/matlab/swz/swz_write_mhm_input.m
new file mode 100644
index 000000000..a687a9494
--- /dev/null
+++ b/matlab/swz/swz_write_mhm_input.m
@@ -0,0 +1,38 @@
+function swz_write_mhm_input(fname,options_ms)
+
+    fh = fopen(fname,'w');
+
+
+    fprintf(fh,'/**********************************************************\n');
+    fprintf(fh,' *** This input file is read by swzmsbvar_mhm_1 and swzmsbvar_mhm_1.exe only, NOT by swzmsbvar_printdraws.exe.\n');
+    fprintf(fh,' ***\n');
+    fprintf(fh,' **********************************************************/\n');
+
+    fprintf(fh,'\n\n//------------- 1st set of posterior draws to find optimal scales for Metropolis (30000). ---------------\n');
+    fprintf(fh,'//== number draws for first burn-in ==//  //For determining the Metropolis scales only.\n');
+    fprintf(fh,'%d\n\n',options_ms.draws_nbr_burn_in_1);
+
+    fprintf(fh,'//------------- MCMC burn-in draws once the Metropolis scales (previous stage) are fixed. --------------\n');
+    fprintf(fh,'//------------- 2nd set of standard burn-in posterior draws to throw away the initial draws (10000). ---------------\n');
+    fprintf(fh,'//== number draws for second burn-in ==//\n');
+    fprintf(fh,'%d\n\n',options_ms.draws_nbr_burn_in_2);
+
+    fprintf(fh,'//--------------- 1st set of posterior draws to compute the mean and variance for the weighting function in the MHM (200000) ----------------\n');
+    fprintf(fh,'//== number draws to estimate mean and variance ==//\n');
+    fprintf(fh,'%d\n\n',options_ms.draws_nbr_mean_var_estimate);
+
+    fprintf(fh,'//--------------- Only applied to mhm_2 process: total number of MCMC draws = thinning factor * 2nd set of saved posterior draws ----------------\n');
+    fprintf(fh,'//== thinning factor for modified harmonic mean process ==//\n');
+    fprintf(fh,'%d\n\n',options_ms.thinning_factor);
+
+    fprintf(fh,'//--------------- 2nd set of saved posterior draws from MHM_2 (second stage): saved draws AFTER thinning (1000000) ----------------\n');
+    fprintf(fh,'//== number draws for modified harmonic mean process ==//\n');
+    fprintf(fh,'%d\n\n',options_ms.draws_nbr_modified_harmonic_mean);
+
+    fprintf(fh,'//------- 1st stage: computing all three tightness factors for Dirichlet.  ---------\n');
+    fprintf(fh,'//------- 2nd stage: hard-code the second scale factor (in principle, we can do all three). ---------\n');
+    fprintf(fh,'//------- It seems that Dan''s code only use the first element of the following scales.  The scale applies to the Dirichlet''s hyperparameter alpha for the diagonal of the transition matrix in the weighting function.  Note that the weighting function for the transition matrix parameters is Dirichlet. ---------\n');
+            
+    fprintf(fh,'//== scale values for Dirichlet distribution ==//\n');
+    fprintf(fh,'1\n\n');
+    fprintf(fh,'%d\n',options_ms.dirichlet_scale);
diff --git a/matlab/swz/sz_prd.m b/matlab/swz/sz_prd.m
index ce05381c9..4ae550001 100644
--- a/matlab/swz/sz_prd.m
+++ b/matlab/swz/sz_prd.m
@@ -1,4 +1,4 @@
-function sz_prd(options_)
+function sz_prd(M,options_)
 %==========================================================================
 %== Directory structure
 %==========================================================================
@@ -53,7 +53,7 @@ markov_file = [options_.ms.markov_file '.dat'];
 %options_.ms.mhm_file = 'MHM_input.dat';
 
 
-mhm_file = [mhm_spec_path '/MHM_input.dat'];
+mhm_file = '/MHM_input.dat';
 %options_.ms.proposal_draws = 100000;
 
 %==========================================================================
@@ -277,7 +277,7 @@ if options_.ms.create_initialization_file == 1
             Ui{kj} = eye(nvar);  Vi{kj} = eye(ncoef);
         end
     else
-        eval(['[Ui,Vi,n0,np,ixmC0Pres] = ' options_.ms.restriction_fname '(options_.ms.nlags,nvar,nexo,indxC0Pres);'])
+        [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
@@ -326,7 +326,7 @@ if options_.ms.create_initialization_file == 1
         %
         tic
         [fhat,xhat,grad,Hhat,itct,fcount,retcodehat] = csminwel('fn_a0freefun',x,H0,'fn_a0freegrad',crit,nit,Ui,nvar,n0,fss,H0inv);
-        endtime = toc
+        endtime = toc;
 
         A0hat = fn_tran_b2a(xhat,Ui,nvar,n0);
 
@@ -480,9 +480,11 @@ if options_.ms.create_initialization_file == 1
     %======================================================================
     %== Create C initialization filename
     %======================================================================
+    swz_write_markov_file(markov_file,M,options_)
     if use_linux == 1
-        create_init_file=[c_path,'/sbvar_init_file ',matlab_filename,' ',m_spec_path,'/',markov_file,' ',options_.ms.output_file_tag];
-        system(create_init_file) %Run operating system command and return result
+        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.exe ',matlab_filename,' ',m_spec_path,'\',markov_file,' ',options_.ms.output_file_tag];
         dos(create_init_file)
@@ -495,10 +497,10 @@ end
 %==========================================================================
 if options_.ms.estimate_msmodel == 1
     if use_linux == 1
-        perform_estimation=[c_path,'/sbvar_estimate -ft ',options_.ms.output_file_tag];
-        system(perform_estimation)
+        perform_estimation=[c_path,'/sbvar_estimation -ft ',options_.ms.output_file_tag];
+        %        system(perform_estimation);
     else
-        perform_estimation=[c_path,'\sbvar_estimate.exe -ft ',options_.ms.output_file_tag];
+        perform_estimation=[c_path,'\sbvar_estimation.exe -ft ',options_.ms.output_file_tag];
         dos(perform_estimation)
     end
 end
@@ -508,13 +510,15 @@ end
 %== Compute marginal data density
 %==========================================================================
 if options_.ms.compute_mdd == 1
+    mhm_file = ['mhm_input_' M.fname '.dat'];
+    swz_write_mhm_input(mhm_file,options_.ms);
     if use_linux == 1
-        compute_mdd1=[c_path,'/sbvar_mhm_1 -ft ',options_.ms.output_file_tag,' -fi ',mhm_spec_path,'/',mhm_file];
+        compute_mdd1=[c_path,'/sbvar_mhm_1 -ft ',options_.ms.output_file_tag,' -fi ',mhm_file];
         system(compute_mdd1);
         compute_mdd2=[c_path,'/sbvar_mhm_2 -ft ',options_.ms.output_file_tag,' -d ',int2str(options_.ms.proposal_draws),' -t 3'];
         system(compute_mdd2);
     else
-        compute_mdd1=[c_path,'\sbvar_mhm_1.exe -ft ',options_.ms.output_file_tag,' -fi ',mhm_spec_path,'/',mhm_file];
+        compute_mdd1=[c_path,'\sbvar_mhm_1.exe -ft ',options_.ms.output_file_tag,' -fi ',mhm_file];
         system(compute_mdd1);
         compute_mdd2=[c_path,'\sbvar_mhm_2.exe -ft ',options_.ms.output_file_tag,' -d ',int2str(options_.ms.proposal_draws),' -t 3'];
         system(compute_mdd2);
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index ff26e0bf2..8cee259a0 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -1100,19 +1100,20 @@ SvarIdentificationStatement::writeOutput(ostream &output, const string &basename
   if (!upper_cholesky_present && !lower_cholesky_present)
     {
       int n = symbol_table.endo_nbr();
-      int m = symbol_table.exo_nbr();
-      int r = getMaxLag();
+//       int m = symbol_table.exo_nbr();
+      int m = 1; // this is the constant, not the shocks
+      int r = getMaxLag()+1;
       int k = r*n+m;
 
       if (k < 1)
-        {
-          cerr << "ERROR: lag = " << r
-               << ", number of endogenous variables = " << n
-               << ", number of exogenous variables = " << m
-               << ". If this is not a logical error in the specification"
-               << " of the .mod file, please report it to the Dynare Team." << endl;
-          exit(EXIT_FAILURE);
-        }
+	{
+	  cerr << "ERROR: lag = " << r
+	       << ", number of endogenous variables = " << n
+	       << ", number of exogenous variables = " << m
+	       << ". If this is not a logical error in the specification"
+	       << " of the .mod file, please report it to the Dynare Team." << endl;
+	  exit(EXIT_FAILURE);
+	}
       if (n < 1)
         {
           cerr << "ERROR: Number of endogenous variables = " << n << "< 1. If this is not a logical "
@@ -1140,7 +1141,7 @@ SvarIdentificationStatement::writeOutput(ostream &output, const string &basename
                 }
 
               if (it->first.first == 0)
-                output << "options_.ms.Qi(" << h+1 << ", " << j << ", "<< i << ");" << endl;
+                output << "options_.ms.Qi(" << h+1 << ", " << j << ", "<< i << ") = 1;" << endl;
               else if (it->first.first > 0)
                 {
                   if ((it->first.first-1)*n+j > k)
@@ -1148,7 +1149,7 @@ SvarIdentificationStatement::writeOutput(ostream &output, const string &basename
                       cerr << "ERROR: lag =" << it->first.first << ", num endog vars = " << n << "current endog var index = " << j << ". Index "
                            << "out of bounds. If the above does not represent a logical error, please report this to the Dyanre Team." << endl;
                     }
-                  output << "options_.ms.Ri(" << h+1 << ", " << (it->first.first-1)*n+j << ", "<< i << ");" << endl;
+                  output << "options_.ms.Ri(" << h+1 << ", " << (it->first.first-1)*n+j << ", "<< i << ") = 1;" << endl;
                 }
               else
                 {
diff --git a/tests/swz/test_upper_cholesky.mod b/tests/swz/test_upper_cholesky.mod
new file mode 100644
index 000000000..74f55ef2a
--- /dev/null
+++ b/tests/swz/test_upper_cholesky.mod
@@ -0,0 +1,16 @@
+addpath '/home/michel/dynare/svn/dynare/trunk/matlab/swz';
+var Y Pie R;
+
+model;
+Y = 0;
+Pie = 0;
+R = 0;
+end;
+
+varobs Y Pie R;
+
+
+sbvar(datafile = data,freq=4,initial_year=1959,final_year=2005,nlags=4,restriction_fname=upper_cholesky);//for SBVAR
+
+//ms_sbvar(datafile = data,freq=4,initial_year=1959,final_year=2005,nlags=4,restriction_fname=ftd_upperchol3v,
+//         markov_file=specification_2v2c,mhm_file=MHM_input); //for markov switching
\ No newline at end of file