diff --git a/matlab/CutSample.m b/matlab/CutSample.m
index c5ba1af60..bfbfeeed9 100644
--- a/matlab/CutSample.m
+++ b/matlab/CutSample.m
@@ -1,10 +1,12 @@
-function CutSample()
+function CutSample(M_, options_, estim_params_)
 
 % function CutSample()
 % Takes a subset from metropolis
 %
 % INPUTS
-%    none
+%   options_         [structure]
+%   estim_params_    [structure]
+%   M_               [structure]
 %
 % OUTPUTS
 %    none
diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m
index eae4aa984..7f6abec49 100644
--- a/matlab/GetPosteriorParametersStatistics.m
+++ b/matlab/GetPosteriorParametersStatistics.m
@@ -1,12 +1,16 @@
-function get_posterior_parameters_statistics()
+function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_)
 % This function prints and saves posterior estimates after the mcmc
 % (+updates of oo_ & TeX output). 
 % 
 % INPUTS 
-%   None.
+%   estim_params_    [structure] 
+%   M_               [structure]
+%   options_         [structure]
+%   bayestopt_       [structure]
+%   oo_              [structure]
 %  
 % OUTPUTS 
-%   None.  
+%   oo_              [structure]  
 %
 % SPECIAL REQUIREMENTS
 %   None.
@@ -30,9 +34,9 @@ function get_posterior_parameters_statistics()
 
 global estim_params_ M_ options_ bayestopt_ oo_
 
-if ~options_.mh_replic & options_.load_mh_file
-   load([M_.fname '_results.mat'],'oo_'); 
-end
+%if ~options_.mh_replic & options_.load_mh_file
+%   load([M_.fname '_results.mat'],'oo_'); 
+%end
 
 TeX   	= options_.TeX;
 nblck 	= options_.mh_nblck;
@@ -60,7 +64,12 @@ tit2 = sprintf('%10s %7s %10s %14s %4s %6s\n',' ','prior mean','post. mean','con
 pformat = '%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f';
 
 disp(' ');disp(' ');disp('ESTIMATION RESULTS');disp(' ');
-disp(sprintf('Log data density is %f.',oo_.MarginalDensity.ModifiedHarmonicMean))
+try
+    disp(sprintf('Log data density is %f.',oo_.MarginalDensity.ModifiedHarmonicMean))
+catch
+    [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_)
+    disp(sprintf('Log data density is %f.',oo_.MarginalDensity.ModifiedHarmonicMean))
+end
 if np
     type = 'parameters';
     if TeX
@@ -78,8 +87,16 @@ if np
             name = bayestopt_.name{ip};
             oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
         else
-            name = bayestopt_.name{ip};
-            [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
+            try
+                name = bayestopt_.name{ip};
+                [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
+            catch
+                Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
+                [post_mean, post_median, post_var, hpd_interval, post_deciles, ...
+                 density] = posterior_moments(Draws,1,options_.mh_conf_sig);
+                name = bayestopt_.name{ip};
+                oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);                
+            end
         end
         disp(sprintf(pformat,name,bayestopt_.pmean(ip),...
                      post_mean, ...
@@ -115,9 +132,19 @@ if nvx
             oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
             M_.Sigma_e(k,k) = post_mean*post_mean;
         else
-            k = estim_params_.var_exo(i,1);
-            name = deblank(M_.exo_names(k,:));
-            [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
+            try
+                k = estim_params_.var_exo(i,1);
+                name = deblank(M_.exo_names(k,:));
+                [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
+            catch
+                Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
+                [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
+                    posterior_moments(Draws,1,options_.mh_conf_sig);
+                k = estim_params_.var_exo(i,1);
+                name = deblank(M_.exo_names(k,:));
+                oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
+                M_.Sigma_e(k,k) = post_mean*post_mean;
+            end
         end
         disp(sprintf(pformat,name,bayestopt_.pmean(ip),post_mean,hpd_interval,...
                      pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
@@ -148,8 +175,16 @@ if nvn
             name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
             oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
         else
-            name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
-            [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
+            try
+                name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
+                [post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
+            catch
+                Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
+                [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
+                    posterior_moments(Draws,1,options_.mh_conf_sig);
+                name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
+                oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
+            end
         end
         disp(sprintf(pformat,name,bayestopt_.pmean(ip),post_mean,hpd_interval, ...
                      pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
@@ -185,11 +220,24 @@ if ncx
             M_.Sigma_e(k1,k2) = post_mean*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
             M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
         else
-            k1 = estim_params_.corrx(i,1);
-            k2 = estim_params_.corrx(i,2);
-            name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
-            NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
-            [post_mean,hpd_interval,post_var] = Extractoo(oo_,NAME,type);
+            try
+                k1 = estim_params_.corrx(i,1);
+                k2 = estim_params_.corrx(i,2);
+                name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
+                NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
+                [post_mean,hpd_interval,post_var] = Extractoo(oo_,NAME,type);
+            catch
+                Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
+                [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
+                    posterior_moments(Draws,1,options_.mh_conf_sig);
+                k1 = estim_params_.corrx(i,1);
+                k2 = estim_params_.corrx(i,2);
+                name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
+                NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
+                oo_ = Filloo(oo_,NAME,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
+                M_.Sigma_e(k1,k2) = post_mean*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
+                M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
+            end
         end
         disp(sprintf(pformat, name,bayestopt_.pmean(ip),post_mean,hpd_interval, ...
                      pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
@@ -224,11 +272,23 @@ if ncn
             oo_ = Filloo(oo_,NAME,type,post_mean,hpd_interval,...
                          post_median,post_var,post_deciles,density);
         else
-            k1 = estim_params_.corrn(i,1);
-            k2 = estim_params_.corrn(i,2);
-            name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];
-            NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
-            [post_mean,hpd_interval,post_var] = Extractoo(oo_,NAME,type);
+            try
+                k1 = estim_params_.corrn(i,1);
+                k2 = estim_params_.corrn(i,2);
+                name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];
+                NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
+                [post_mean,hpd_interval,post_var] = Extractoo(oo_,NAME,type);
+            catch
+                Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
+                [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
+                    posterior_moments(Draws,1,options_.mh_conf_sig);
+                k1 = estim_params_.corrn(i,1);
+                k2 = estim_params_.corrn(i,2);
+                name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];
+                NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
+                oo_ = Filloo(oo_,NAME,type,post_mean,hpd_interval,...
+                             post_median,post_var,post_deciles,density);
+            end
         end
         disp(sprintf(pformat, name,bayestopt_.pmean(ip),post_mean,hpd_interval, ...
                      pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
diff --git a/matlab/McMCDiagnostics.m b/matlab/McMCDiagnostics.m
index aa3d56a5c..98c0923cb 100644
--- a/matlab/McMCDiagnostics.m
+++ b/matlab/McMCDiagnostics.m
@@ -1,10 +1,12 @@
-function McMCDiagnostics
+function McMCDiagnostics(options_, estim_params_, M_)
 % function McMCDiagnostics
 % Computes convergence tests 
 % 
 % INPUTS 
-%   none
-%  
+%   options_         [structure]
+%   estim_params_    [structure]
+%   M_               [structure]
+%
 % OUTPUTS 
 %   none  
 %
@@ -28,8 +30,6 @@ function McMCDiagnostics
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global options_ estim_params_ M_ 
-
 DirectoryName = CheckPath('Output');
 MhDirectoryName = CheckPath('metropolis');
 
diff --git a/matlab/PlotPosteriorDistributions.m b/matlab/PlotPosteriorDistributions.m
index a7f68ad9d..1d1cac5e5 100644
--- a/matlab/PlotPosteriorDistributions.m
+++ b/matlab/PlotPosteriorDistributions.m
@@ -1,13 +1,17 @@
-function PlotPosteriorDistributions()
+function oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt_, oo_)
 
 % function PlotPosteriorDistributions()
 % plots posterior distributions
 %
 % INPUTS
-%    none
+%    estim_params_   [structure] 
+%    M_              [structure]
+%    options_        [structure] 
+%    bayestopt_      [structure]
+%    oo_             [structure]
 %    
 % OUTPUTS
-%    none
+%    oo_             [structure]  
 %    
 % SPECIAL REQUIREMENTS
 %    none
diff --git a/matlab/compute_moments_varendo.m b/matlab/compute_moments_varendo.m
new file mode 100644
index 000000000..3b68842d5
--- /dev/null
+++ b/matlab/compute_moments_varendo.m
@@ -0,0 +1,57 @@
+function oo_ = compute_moments_varendo(options_,M_,oo_,var_list_)
+% Computes the second order moments (autocorrelation function, covariance
+% matrix and variance decomposition) for all the endogenous variables selected in
+% var_list_. The results are saved in oo_
+%  
+% INPUTS:
+%   options_        [structure]    Dynare structure.
+%   M_              [structure]    Dynare structure (related to model definition).
+%   oo_             [structure]    Dynare structure (results).
+%   var_list_       [string]       Array of string with endogenous variable names.
+%    
+% OUTPUTS
+%   oo_             [structure]    Dynare structure (results).
+%
+% SPECIAL REQUIREMENTS
+%   none
+
+% Copyright (C) 2008 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/>.    
+    NumberOfEndogenousVariables = rows(var_list_);
+    NumberOfExogenousVariables = M_.exo_nbr;
+    list_of_exogenous_variables = M_.exo_names;
+    NumberOfLags = options_.ar;
+    % COVARIANCE MATRIX.
+    for i=1:NumberOfEndogenousVariables
+        for j=i:NumberOfEndogenousVariables
+            oo_ = posterior_analysis('variance',var_list_(i,:),var_list_(j,:),[],options_,M_,oo_);
+        end
+    end
+    % CORRELATION FUNCTION.
+    for h=NumberOfLags:-1:1
+        for i=1:NumberOfEndogenousVariables
+            for j=1:NumberOfEndogenousVariables
+                oo_ = posterior_analysis('correlation',var_list_(i,:),var_list_(j,:),h,options_,M_,oo_);
+            end
+        end
+    end
+    % VARIANCE DECOMPOSITION.
+    for i=1:NumberOfEndogenousVariables
+        for j=i:NumberOfExogenousVariables
+            oo_ = posterior_analysis('decomposition',var_list_(i,:),M_.exo_names(j,:),[],options_,M_,oo_);
+        end
+    end
\ No newline at end of file
diff --git a/matlab/correlation_posterior_analysis.m b/matlab/correlation_posterior_analysis.m
index f7d21705a..f7cefec01 100644
--- a/matlab/correlation_posterior_analysis.m
+++ b/matlab/correlation_posterior_analysis.m
@@ -76,23 +76,38 @@ function oo_ = correlation_posterior_analysis(SampleSize,dname,fname,vartan,nvar
         tmp(i1:i2) = Correlation_array(:,indx1,indx2,nar);
         i1 = i2+1;
     end
-    [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
-        posterior_moments(tmp,1,mh_conf_sig);
     name = [ var1 '.' var2 ];
-    if isfield(oo_,'PosteriorTheoreticalMoments')
-        if isfield(oo_.PosteriorTheoreticalMoments,'dsge')
-            if isfield(oo_.PosteriorTheoreticalMoments.dsge,'correlation')
-                oo_ = fill_output_structure(var1,var2,oo_,'mean',nar,post_mean);
-                oo_ = fill_output_structure(var1,var2,oo_,'median',nar,post_median);
-                oo_ = fill_output_structure(var1,var2,oo_,'variance',nar,post_var);
-                oo_ = fill_output_structure(var1,var2,oo_,'hpdinf',nar,hpd_interval(1));
-                oo_ = fill_output_structure(var1,var2,oo_,'hpdsup',nar,hpd_interval(2));
-                oo_ = fill_output_structure(var1,var2,oo_,'deciles',nar,post_deciles);
-                oo_ = fill_output_structure(var1,var2,oo_,'density',nar,density);
+    if ~isconst(tmp)
+        [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
+            posterior_moments(tmp,1,mh_conf_sig);
+        if isfield(oo_,'PosteriorTheoreticalMoments')
+            if isfield(oo_.PosteriorTheoreticalMoments,'dsge')
+                if isfield(oo_.PosteriorTheoreticalMoments.dsge,'correlation')
+                    oo_ = fill_output_structure(var1,var2,oo_,'mean',nar,post_mean);
+                    oo_ = fill_output_structure(var1,var2,oo_,'median',nar,post_median);
+                    oo_ = fill_output_structure(var1,var2,oo_,'variance',nar,post_var);
+                    oo_ = fill_output_structure(var1,var2,oo_,'hpdinf',nar,hpd_interval(1));
+                    oo_ = fill_output_structure(var1,var2,oo_,'hpdsup',nar,hpd_interval(2));
+                    oo_ = fill_output_structure(var1,var2,oo_,'deciles',nar,post_deciles);
+                    oo_ = fill_output_structure(var1,var2,oo_,'density',nar,density);
+                end
+            end
+        end
+    else
+        if isfield(oo_,'PosteriorTheoreticalMoments')
+            if isfield(oo_.PosteriorTheoreticalMoments,'dsge')
+                if isfield(oo_.PosteriorTheoreticalMoments.dsge,'correlation')
+                    oo_ = fill_output_structure(var1,var2,oo_,'mean',nar,NaN);
+                    oo_ = fill_output_structure(var1,var2,oo_,'median',nar,NaN);
+                    oo_ = fill_output_structure(var1,var2,oo_,'variance',nar,NaN);
+                    oo_ = fill_output_structure(var1,var2,oo_,'hpdinf',nar,NaN);
+                    oo_ = fill_output_structure(var1,var2,oo_,'hpdsup',nar,NaN);
+                    oo_ = fill_output_structure(var1,var2,oo_,'deciles',nar,NaN);
+                    oo_ = fill_output_structure(var1,var2,oo_,'density',nar,NaN);
+                end
             end
         end
     end
-
     
 function oo_ = initialize_output_structure(var1,var2,nar,oo_)
     name = [ var1 '.' var2 ];
diff --git a/matlab/covariance_posterior_analysis.m b/matlab/covariance_posterior_analysis.m
index 98921eaca..b879c5731 100644
--- a/matlab/covariance_posterior_analysis.m
+++ b/matlab/covariance_posterior_analysis.m
@@ -49,12 +49,12 @@ function oo_ = covariance_posterior_analysis(NumberOfSimulations,dname,fname,var
                             return
                         end
                     end
-                end         
+                end
             end
         end
     end
     tmp = dir([ dname '/metropolis/'  fname '_Posterior2ndOrderMoments*.mat']);
-    NumberOfFiles= length(tmp);
+    NumberOfFiles = length(tmp);
     i1 = 1; tmp = zeros(NumberOfSimulations,1);
     for file = 1:NumberOfFiles
         load([ dname '/metropolis/'  fname '_Posterior2ndOrderMoments' int2str(file) '.mat']);
@@ -62,13 +62,23 @@ function oo_ = covariance_posterior_analysis(NumberOfSimulations,dname,fname,var
         tmp(i1:i2) = Covariance_matrix(:,symmetric_matrix_index(indx1,indx2,nvar));
         i1 = i2+1;
     end
-    [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
-        posterior_moments(tmp,1,mh_conf_sig);
     name = [var1 '.' var2];
-    eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.mean.' name ' = post_mean;']);
-    eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.median.' name ' = post_median;']);
-    eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.variance.' name ' = post_var;']);
-    eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.hpdinf.' name ' = hpd_interval(1);']);
-    eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.hpdsup.' name ' = hpd_interval(2);']);
-    eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.deciles.' name ' = post_deciles;']);
-    eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.density.' name ' = density;']);
\ No newline at end of file
+    if ~isconst(tmp)
+        [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
+            posterior_moments(tmp,1,mh_conf_sig);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.mean.' name ' = post_mean;']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.median.' name ' = post_median;']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.variance.' name ' = post_var;']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.hpdinf.' name ' = hpd_interval(1);']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.hpdsup.' name ' = hpd_interval(2);']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.deciles.' name ' = post_deciles;']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.density.' name ' = density;']);
+    else
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.mean.' name ' = NaN;']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.median.' name ' = NaN;']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.variance.' name ' = NaN;']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.hpdinf.' name ' = NaN;']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.hpdsup.' name ' = NaN;']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.deciles.' name ' = NaN;']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.density.' name ' = NaN;']);
+    end
\ No newline at end of file
diff --git a/matlab/dsge_posterior_theoretical_covariance.m b/matlab/dsge_posterior_theoretical_covariance.m
index 3b238cf8d..cb4e58c50 100644
--- a/matlab/dsge_posterior_theoretical_covariance.m
+++ b/matlab/dsge_posterior_theoretical_covariance.m
@@ -40,6 +40,8 @@ type = 'posterior';
     
 % Set varlist (vartan)
 [ivar,vartan] = set_stationary_variables_list;
+vartan
+pause
 nvar = length(ivar);
 
 % Set the size of the auto-correlation function to zero.
@@ -68,7 +70,7 @@ CovarFileNumber = 1;
 % Compute 2nd order moments and save them in *_Posterior2ndOrderMoments* files
 linea = 0;
 for file = 1:NumberOfDrawsFiles
-    load([M_.dname '/metropolis/' DrawsFiles(file).name '.mat']);
+    load([M_.dname '/metropolis/' DrawsFiles(file).name ],'pdraws');
     NumberOfDraws = rows(pdraws);
     isdrsaved = cols(pdraws)-1;
     for linee = 1:NumberOfDraws
diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m
index b6b638aef..4be555372 100644
--- a/matlab/dynare_estimation.m
+++ b/matlab/dynare_estimation.m
@@ -616,7 +616,7 @@ if any(bayestopt_.pshape > 0) & options_.posterior_mode_estimation
       eval(['oo_.posterior_std.measurement_errors_corr.' NAME ' = stdh(ip);']); 
       ip = ip+1;
     end
-  end  
+  end
   %% Laplace approximation to the marginal log density:
   if ~options_.bvar_dsge
     md_Laplace = .5*size(xparam1,1)*log(2*pi) + .5*log(det(invhess)) ...
@@ -898,16 +898,16 @@ end
 
 if (any(bayestopt_.pshape  >0 ) & options_.mh_replic) | ...
       (any(bayestopt_.pshape >0 ) & options_.load_mh_file)  %% not ML estimation
-  bounds = prior_bounds(bayestopt_);
-  bounds(:,1)=max(bounds(:,1),lb);
-  bounds(:,2)=min(bounds(:,2),ub);
-  bayestopt_.lb = bounds(:,1);
-  bayestopt_.ub = bounds(:,2);
-  if any(xparam1 < bounds(:,1)) | any(xparam1 > bounds(:,2))
-    find(xparam1 < bounds(:,1))
-    find(xparam1 > bounds(:,2))
-    error('Mode values are outside prior bounds. Reduce prior_trunc.')
-  end
+    bounds = prior_bounds(bayestopt_);
+    bounds(:,1)=max(bounds(:,1),lb);
+    bounds(:,2)=min(bounds(:,2),ub);
+    bayestopt_.lb = bounds(:,1);
+    bayestopt_.ub = bounds(:,2);
+    if any(xparam1 < bounds(:,1)) | any(xparam1 > bounds(:,2))
+        find(xparam1 < bounds(:,1))
+        find(xparam1 > bounds(:,2))
+        error('Mode values are outside prior bounds. Reduce prior_trunc.')
+    end
   % runs MCMC
   if options_.mh_replic
     if options_.load_mh_file & options_.use_mh_covariance_matrix
@@ -920,26 +920,29 @@ if (any(bayestopt_.pshape  >0 ) & options_.mh_replic) | ...
     end
   end
   if ~options_.nodiagnostic & options_.mh_replic > 1000 & options_.mh_nblck > 1
-    McMCDiagnostics;
+      McMCDiagnostics(options_, estim_params_, M_);
   end
   %% Here i discard first half of the draws:
   CutSample;
   %% Estimation of the marginal density from the Mh draws:
   if options_.mh_replic
-      marginal = marginal_density;
+      [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_);
   end
   %%
-  GetPosteriorParametersStatistics;
+  oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_);
   %% Results are saved (in case of an anormal exit from dynare or matlab)...
-  save([M_.fname '_results.mat'],'oo_','M_');
+  %%save([M_.fname '_results.mat'],'oo_','M_');
   %%
-  PlotPosteriorDistributions;
+  oo_
+  oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt_, oo_);
   metropolis_draw(1);
   if options_.bayesian_irf
     PosteriorIRF('posterior');
   end
-  if options_.smoother | ~isempty(options_.filter_step_ahead) | options_.forecast ...
-          | options_.moments_varendo
+  if options_.moments_varendo
+      oo_ = compute_moments_varendo(options_,M_,oo_,var_list_);
+  end
+  if options_.smoother | ~isempty(options_.filter_step_ahead) | options_.forecast
     prior_posterior_statistics('posterior',data,gend);
   end
 
diff --git a/matlab/isconst.m b/matlab/isconst.m
new file mode 100644
index 000000000..9bcea59bd
--- /dev/null
+++ b/matlab/isconst.m
@@ -0,0 +1,32 @@
+function aa = isconst(y)
+% Returns 1 if vector y is constant, 0 otherwise.
+%  
+% INPUTS:
+%   yy        [double]    n*1 vector.
+%
+% OUTPUTS
+%   aa        [integer]   scalar equal to 1 or 0.
+%
+% SPECIAL REQUIREMENTS
+%   none
+
+% Copyright (C) 2008 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/>
+    aa = 0;
+    if all(abs(y(2:end)-y(1:end-1))<1e-10)
+        aa = 1;
+    end
\ No newline at end of file
diff --git a/matlab/marginal_density.m b/matlab/marginal_density.m
index 17481681d..a7074b997 100644
--- a/matlab/marginal_density.m
+++ b/matlab/marginal_density.m
@@ -1,12 +1,16 @@
-function marginal = marginal_density()
+function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_)
 % function marginal = marginal_density()
 % Computes the marginal density
 %
 % INPUTS
-%    none
+%   options_         [structure]
+%   estim_params_    [structure]
+%   M_               [structure]
+%   oo_              [structure]
 %
 % OUTPUTS
-%    marginal:  marginal density
+%   marginal:        [double]     marginal density (modified harmonic mean)
+%   oo_              [structure]
 %
 % SPECIAL REQUIREMENTS
 %    none
@@ -28,7 +32,6 @@ function marginal = marginal_density()
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global M_ options_ estim_params_ oo_
 
 npar = estim_params_.np+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nvx;
 nblck = options_.mh_nblck;
diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m
index fd7efea84..a0f522cad 100644
--- a/matlab/prior_posterior_statistics.m
+++ b/matlab/prior_posterior_statistics.m
@@ -50,7 +50,7 @@ exo_nbr=M_.exo_nbr;
 nvobs     = size(options_.varobs,1);
 iendo = 1:endo_nbr;
 horizon = options_.forecast;
-moments_varendo = options_.moments_varendo;
+% moments_varendo = options_.moments_varendo;
 filtered_vars = options_.filtered_vars;
 if horizon
     i_last_obs = gend+(1-M_.maximum_endo_lag:0);
@@ -107,8 +107,8 @@ else
   end
 end
 
-irun = ones(8,1);
-ifil = zeros(8,1);
+irun = ones(7,1);
+ifil = zeros(7,1);
 
 if exist('OCTAVE_VERSION')
     diary off;
@@ -138,16 +138,16 @@ if options_.forecast
     stock_forcst_point = zeros(endo_nbr,horizon+maxlag,MAX_nforc2);
     run_smoother = 1;
 end
-if moments_varendo
-    stock_moments = cell(MAX_momentsno,1);
-end
+%if moments_varendo
+%    stock_moments = cell(MAX_momentsno,1);
+%end
 for b=1:B
   [deep, logpo] = GetOneDraw(type);
   set_all_parameters(deep);
   dr = resol(oo_.steady_state,0);
-  if moments_varendo
-      stock_moments{irun(8)} = compute_model_moments(dr,M_,options_);
-  end
+  %if moments_varendo
+  %    stock_moments{irun(8)} = compute_model_moments(dr,M_,options_);
+  %end
   if run_smoother
       [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK] = ...
           DsgeSmoother(deep,gend,Y);
@@ -205,7 +205,7 @@ for b=1:B
   stock_logpo(irun(5),1) = logpo;
   stock_ys(irun(5),:) = SteadyState';
 
-  irun = irun +  ones(8,1);
+  irun = irun +  ones(7,1);
 
   if irun(1) > MAX_nsmoo || b == B
       stock = stock_smooth(:,:,1:irun(1)-1);
@@ -258,12 +258,12 @@ for b=1:B
       irun(7) = 1;
   end
 
-  if moments_varendo && (irun(8) > MAX_momentsno || b == B)
-      stock = stock_moments(1:irun(8)-1);
-      ifil(8) = ifil(8) + 1;
-      save([DirectoryName '/' M_.fname '_moments' int2str(ifil(8)) '.mat'],'stock');
-      irun(8) = 1;
-  end
+  % if moments_varendo && (irun(8) > MAX_momentsno || b == B)
+  %    stock = stock_moments(1:irun(8)-1);
+  %    ifil(8) = ifil(8) + 1;
+  %    save([DirectoryName '/' M_.fname '_moments' int2str(ifil(8)) '.mat'],'stock');
+  %    irun(8) = 1;
+  % end
   
   if exist('OCTAVE_VERSION')
       printf('Taking subdraws: %3.f%% done\r', b/B);