diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 1b2fdc19b3b000e4d757b16862daf4aa90aeb000..65a5d4ec3aba3b051410a01e3b72c1b27b7fe4ae 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -339,7 +339,8 @@ if ~options_.mh_posterior_mode_estimation
         oo_.posterior.optimization.log_density=-fval;
     end
     if options_.cova_compute
-        invhess = inv(hh);
+        hsd = sqrt(diag(hh));
+        invhess = inv(hh./(hsd*hsd'))./(hsd*hsd');
         stdh = sqrt(diag(invhess));
         oo_.posterior.optimization.Variance = invhess;
     end
@@ -365,8 +366,7 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
     % Laplace approximation to the marginal log density:
     if options_.cova_compute
         estim_params_nbr = size(xparam1,1);
-        scale_factor = -sum(log10(diag(invhess)));
-        log_det_invhess = -estim_params_nbr*log(scale_factor)+log(det(scale_factor*invhess));
+        log_det_invhess = log(det(invhess./(stdh*stdh')))+2*sum(log(stdh));
         likelihood = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
         oo_.MarginalDensity.LaplaceApproximation = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess - likelihood;
         skipline()
diff --git a/tests/estimation/fs2000.mod b/tests/estimation/fs2000.mod
index 640f1961f96b49a1d5f6fc79ae497e3fd31ade69..a4d4030b2ca2a35113ddb136b0809847ed8aae38 100644
--- a/tests/estimation/fs2000.mod
+++ b/tests/estimation/fs2000.mod
@@ -91,8 +91,9 @@ copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat'],[M_.d
 estimation(mode_compute=0,mode_file=fs2000_mode,order=1, datafile=fsdat_simul, nobs=192, loglinear, mh_replic=1500, mh_nblocks=1, mh_jscale=0.8);
 hh=eye(size(bayestopt_.name,1));
 save('fs2000_mode','hh','-append')
+Laplace = oo_.MarginalDensity.LaplaceApproximation;
 estimation(mode_compute=0,mode_file=fs2000_mode,order=1, datafile=fsdat_simul, nobs=192, loglinear, mh_replic=1500, mh_nblocks=1, mh_jscale=10,load_mh_file);
-
+oo_.MarginalDensity.LaplaceApproximation = Laplace;
 
 temp1=load([M_.dname '_mh1_blck1.mat']);
 temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat']);
diff --git a/tests/estimation/fs2000_model_comparison.mod b/tests/estimation/fs2000_model_comparison.mod
index 7f2a3253a1d330d63a5a50499b42d99b9161d51d..c88135d16995c4f33cf01d732965c63f56684f2f 100644
--- a/tests/estimation/fs2000_model_comparison.mod
+++ b/tests/estimation/fs2000_model_comparison.mod
@@ -99,7 +99,7 @@ if oo_.Model_Comparison.fs2000.Posterior_Model_Probability < oo_.Model_Compariso
 end
 oo_laplace=oo_;
 model_comparison (marginal_density=modifiedharmonicmean) fs2000(0.5) fs2000_initialize_from_calib(0.75);
-if abs(oo_laplace.Model_Comparison.fs2000.Log_Marginal_Density-oo_.Model_Comparison.fs2000.Log_Marginal_Density)>0.2
+if abs(oo_laplace.Model_Comparison.fs2000.Log_Marginal_Density-oo_.Model_Comparison.fs2000.Log_Marginal_Density)>0.5
    error('Laplace and Harmonic Mean do not match') 
 end
 model_comparison (marginal_density=modifiedharmonicmean) fs2000(0) fs2000_initialize_from_calib(1);