diff --git a/matlab/compute_mh_covariance_matrix.m b/matlab/compute_mh_covariance_matrix.m
index a11050909..ed809035b 100644
--- a/matlab/compute_mh_covariance_matrix.m
+++ b/matlab/compute_mh_covariance_matrix.m
@@ -1,12 +1,13 @@
-function covar = compute_mh_covariance_matrix()
-% Estimation of the posterior covariance matrix. 
+function [m0,s0] = compute_mh_covariance_matrix()
+% Estimation of the posterior covariance matrix and expectation. 
 % 
 % INPUTS 
 %   None.
 %  
-% OUTPUTS 
-%   o covar   [double] p*p matrix, posterior covariance of the estimated
-%   parameters (computed from previous metropolis hastings).  
+% OUTPUTS
+%   o  m0  [double]  (n*1) vector, posterior expectation of the parameters.
+%   o  s0  [double]  (n*n) matrix, posterior covariance of the parameters 
+%                    (computed from previous metropolis hastings).
 %
 %
 % ALGORITHM 
@@ -18,43 +19,49 @@ function covar = compute_mh_covariance_matrix()
 %  
 % part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
 % Gnu Public License.
-global M_ options_ estim_params_ oo_
+global M_ options_ estim_params_
 
-npar = estim_params_.np+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.nvx;
+n = 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;
+FirstLine = record.KeepedDraws.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);
+load([ MhDirectoryName '/' M_.fname '_mh' int2str(1) '_blck' int2str(1)],'x2','logpo2');
+params = x2(1,:); oldlogpo2 = logpo2(1);
+for blck = 2:nblck
+    load([ MhDirectoryName '/' M_.fname '_mh' int2str(1) '_blck' int2str(blck)],'x2','logpo2');
+    if logpo2(1)>oldlogpo2
+        oldlogpo2 = logpo2(1);
+        params = x2(1,:);
+    end
+end
 
-fprintf('MH: I''m computing the posterior mean... ');
+offset = 0;
+m0 = zeros(n,1); s0 = zeros(n,n);
 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,:));
+    [tmp,idx] = max(logpo2);
+    if tmp>oldlogpo2
+        oldlogpo2 = tmp;
+        params = x2(idx,:);
+    end
+    [m0,s0,offset] = recursive_moments(m0,s0,x2(FirstLine,:),offset);
   end
-  ifil = 1;
+  FirstLine = 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!)
+
+xparam1 = params';
+hh = inv(s0);
+fval = oldlogpo2;
+
+save([M_.fname 'mh_mode'],'xparam1','hh','fval');
\ No newline at end of file