diff --git a/matlab/PlotPosteriorDistributions.m b/matlab/PlotPosteriorDistributions.m
index 93c98455f275096bda51da3cf8bb8644b49b3d3a..aefa2a5d1292195a32fc6754240c01201dc5ea9f 100644
--- a/matlab/PlotPosteriorDistributions.m
+++ b/matlab/PlotPosteriorDistributions.m
@@ -55,12 +55,16 @@ for i=1:npar
     name = deblank(M_.exo_names(estim_params_.var_exo(i,1),:));  
     eval(['x1 = oo_.posterior_density.shocks_std.' name '(:,1);'])
     eval(['f1 = oo_.posterior_density.shocks_std.' name '(:,2);'])
-    eval(['pmode = oo_.posterior_mode.shocks_std.' name ';'])
+    if options_.posterior_mode_estimation
+      eval(['pmode = oo_.posterior_mode.shocks_std.' name ';'])
+    end
   elseif i <= nvx+nvn
     name = deblank(options_.varobs(estim_params_.var_endo(i-nvx,1),:));
     eval(['x1 = oo_.posterior_density.measurement_errors_std.' name '(:,1);'])
     eval(['f1 = oo_.posterior_density.measurement_errors_std.' name '(:,2);'])    
-    eval(['pmode = oo_.posterior_mode.measurement_errors_std.' name ';'])  
+    if options_.posterior_mode_estimation
+      eval(['pmode = oo_.posterior_mode.measurement_errors_std.' name ';'])
+    end     
   elseif i <= nvx+nvn+ncx
     j = i - (nvx+nvn)
     k1 = estim_params_.corrx(j,1);
@@ -68,7 +72,9 @@ for i=1:npar
     name = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];  
     eval(['x1 = oo_.posterior_density.shocks_corr.' name '(:,1);'])
     eval(['f1 = oo_.posterior_density.shocks_corr.' name '(:,2);'])    
-    eval(['pmode = oo_.posterior_mode.shocks_corr.' name ';'])  
+    if options_.posterior_mode_estimation
+      eval(['pmode = oo_.posterior_mode.shocks_corr.' name ';'])  
+    end
   elseif i <= nvx+nvn+ncx+ncn
     j = i - (nvx+nvn+ncx);
     k1 = estim_params_.corrn(j,1);
@@ -76,13 +82,17 @@ for i=1:npar
     name = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
     eval(['x1 = oo_.posterior_density.measurement_errors_corr.' name '(:,1);'])
     eval(['f1 = oo_.posterior_density.measurement_errors_corr.' name '(:,2);'])
-    eval(['pmode = oo_.posterior_mode.measurement_errors_corr.' name ';'])
+    if options_.posterior_mode_estimation
+      eval(['pmode = oo_.posterior_mode.measurement_errors_corr.' name ';'])
+    end
   else
     j = i - (nvx+nvn+ncx+ncn);
     name = deblank(M_.param_names(estim_params_.param_vals(j,1),:));
     eval(['x1 = oo_.posterior_density.' name '(:,1);'])
     eval(['f1 = oo_.posterior_density.' name '(:,2);'])
-    eval(['pmode = oo_.posterior_mode.parameters.' name ';'])
+    if options_.posterior_mode_estimation
+      eval(['pmode = oo_.posterior_mode.parameters.' name ';'])
+    end
   end
   top1 = max(f1);
   top0 = max([top1;top2]);
@@ -95,7 +105,9 @@ for i=1:npar
   set(hh,'color',[0.7 0.7 0.7]);
   hold on;
   plot(x1,f1,'-k','linewidth',2);
-  plot( [pmode pmode], [0.0 1.1*top0], '--g', 'linewidth', 2);
+  if options_.posterior_mode_estimation
+    plot( [pmode pmode], [0.0 1.1*top0], '--g', 'linewidth', 2);
+  end
   box on;
   axis([borneinf bornesup 0 1.1*top0]);
   title(nam,'Interpreter','none');
diff --git a/matlab/compute_mh_covariance_matrix.m b/matlab/compute_mh_covariance_matrix.m
new file mode 100644
index 0000000000000000000000000000000000000000..a11050909b57f96f1b188668012e79113f123be5
--- /dev/null
+++ b/matlab/compute_mh_covariance_matrix.m
@@ -0,0 +1,60 @@
+function covar = compute_mh_covariance_matrix()
+% Estimation of the posterior covariance matrix. 
+% 
+% INPUTS 
+%   None.
+%  
+% OUTPUTS 
+%   o covar   [double] p*p matrix, posterior covariance of the estimated
+%   parameters (computed from previous metropolis hastings).  
+%
+%
+% ALGORITHM 
+%   None.       
+%
+% SPECIAL REQUIREMENTS
+%   None.
+%  
+%  
+% part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
+% Gnu Public License.
+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;
+
+MhDirectoryName = CheckPath('metropolis');
+load([ MhDirectoryName '/'  M_.fname '_mh_history'])
+
+FirstMhFile = record.KeepedDraws.FirstMhFile;
+FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine;
+TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
+TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
+MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
+TODROP = floor(options_.mh_drop*TotalNumberOfMhDraws);
+
+MU = zeros(1,npar);
+SIGMA = zeros(npar,npar);
+
+fprintf('MH: I''m computing the posterior mean... ');
+for n = FirstMhFile:TotalNumberOfMhFiles
+  for b = 1:nblck
+    load([ MhDirectoryName '/' M_.fname '_mh' int2str(n) '_blck' int2str(b)],'x2','logpo2'); 
+    MU = MU + sum(x2(ifil:end,:));
+  end
+  ifil = 1;
+end
+MU = MU/((TotalNumberOfMhDraws-TODROP)*nblck);
+MU1 = repmat(MU,MAX_nruns,1);	
+fprintf(' Done!\n');
+fprintf('MH: I''m computing the posterior covariance matrix... ');
+ifil = FirstLine;
+for n = FirstMhFile:TotalNumberOfMhFiles
+  for b = 1:nblck
+    load([MhDirectoryName '/' M_.fname '_mh' int2str(n) '_blck' int2str(b)],'x2');
+    x = x2(ifil:end,:)-MU1(1:size(x2(ifil:end,:),1),:);
+    SIGMA = SIGMA + x'*x;
+  end
+  ifil = 1;
+end
+SIGMA =  SIGMA/((TotalNumberOfMhDraws-TODROP)*nblck);%<=== Variance of the parameters (ok!)
diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m
index b42bb8d0b62fa45f7d75d4e2f824e674b5c64c5d..43afba1b871413d381c1d4c5a39f16ac96e4f7f2 100644
--- a/matlab/dynare_estimation.m
+++ b/matlab/dynare_estimation.m
@@ -66,6 +66,7 @@ options_ = set_default_option(options_,'Opt6Iter',3);
 options_ = set_default_option(options_,'Opt6Numb',100000);
 options_ = set_default_option(options_,'steadystate_flag',0);
 options_ = set_default_option(options_,'logdata',0);
+options_ = set_default_option(options_,'use_mh_covariance_matrix',0);
 
 if exist([M_.fname '_steadystate.m'])
   options_.steadystate_flag = 1;
@@ -315,8 +316,11 @@ if options_.mode_compute > 0 & options_.posterior_mode_estimation
         % covariance (a diagonal matrix) % Except for infinite variances ;-)
       stdev = bayestopt_.pstdev;
       indx = find(isinf(stdev));
-      stdev(indx) = sqrt(2)*ones(length(indx),1)*0.1;
+      stdev(indx) = ones(length(indx),1)*0.1;
+      indx = find(stdev>2);
+      stdev(indx) = ones(length(indx),1)*0.1;      
       CovJump = diag(stdev).^2;
+      CovJump = eye(length(stdev))*0.5;
     end
     OldPostVar = CovJump;
     Scale = options_.mh_jscale;
@@ -379,7 +383,9 @@ if options_.posterior_mode_estimation
   invhess = inv(hh);
   stdh = sqrt(diag(invhess));
 else
-  invhess = eye(length(xparam1));
+  variances = bayestopt_.pstdev.^2;
+  invhess = 0.001*diag(variances);
+  invhess = 0.001*eye(length(variances));
 end
   
 if any(bayestopt_.pshape > 0) & options_.posterior_mode_estimation
@@ -756,7 +762,14 @@ if (any(bayestopt_.pshape  >0 ) & options_.mh_replic) | ...
     error('Mode values are outside prior bounds. Reduce prior_trunc.')
   end  
   if options_.mh_replic
-    metropolis('DsgeLikelihood',xparam1,invhess,bounds,gend,data);
+    if ~options_.load_mh_file
+      metropolis('DsgeLikelihood',xparam1,invhess,bounds,gend,data);
+    else
+      if options_.use_mh_covariance_matrix
+        invhess = compute_mh_covariance_matrix();
+      end
+      metropolis('DsgeLikelihood',xparam1,invhess,bounds,gend,data);
+    end
   end
   if ~options_.nodiagnostic & options_.mh_replic > 1000 & options_.mh_nblck > 1
     McMCDiagnostics;
@@ -773,7 +786,6 @@ if (any(bayestopt_.pshape  >0 ) & options_.mh_replic) | ...
     PosteriorIRF('posterior');
   end
   return
-  
 end
 
 if ~((any(bayestopt_.pshape > 0) & options_.mh_replic) | (any(bayestopt_.pshape ...
diff --git a/matlab/gmhmaxlik.m b/matlab/gmhmaxlik.m
index 9af487814d4cc214d1049fe02626f175aeab1c5d..a9724c16655686623ef9e98b7eacf6e9da92e4e2 100644
--- a/matlab/gmhmaxlik.m
+++ b/matlab/gmhmaxlik.m
@@ -208,6 +208,8 @@ if strcmpi(info,'LastCall')
   %%
   %% Now I climb the hill
   %%
+  climb = 0;
+  if climb
   hh = waitbar(0,' ');
   set(hh,'Name','Now I am climbing the hill...')
   j = 1; jj  = 1;
@@ -247,6 +249,7 @@ if strcmpi(info,'LastCall')
     jj = jj + 1;
   end
   close(hh);
+  end%climb
 else
   Scale = iScale;
 end