diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index 5bd3ea664d9174275500897187172231574806ac..b8bd7d89e363fbccc67fc757b03ebc08a588b18a 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -5848,9 +5848,11 @@ block decomposition of the model (see :opt:`block`).
        crashed chain with the respective last random number generator
        state is currently not supported.
 
-    .. option:: mh_mode = INTEGER
+    .. option:: mh_posterior_mode_estimation
 
-        ...
+       Skip optimizer-based mode-finding and instead compute the mode based 
+       on a run of a MCMC. The MCMC will start at the prior mode and use the prior
+       variances to compute the inverse Hessian.
 
     .. option:: mode_file = FILENAME
 
diff --git a/matlab/check_posterior_sampler_options.m b/matlab/check_posterior_sampler_options.m
index 0cfff350d2dd38501df1bc8a4915807d67ae5128..2f7a80b828786cab847b383f73e961b998323ef3 100644
--- a/matlab/check_posterior_sampler_options.m
+++ b/matlab/check_posterior_sampler_options.m
@@ -308,7 +308,7 @@ if init
         % slice posterior sampler does not require mode or hessian to run
         % needs to be set to 1 to skip parts in dynare_estimation_1.m
         % requiring posterior maximization/calibrated smoother before MCMC
-        options_.mh_posterior_mode_estimation=1;
+        options_.mh_posterior_mode_estimation=true;
 
         if ~ posterior_sampler_options.slice_initialize_with_mode
             % by default, slice sampler should trigger
@@ -324,7 +324,7 @@ if init
                 disp('check_posterior_sampler_options:: but no mode file nor posterior maximization is selected,')
                 error('check_posterior_sampler_options:: The option "slice_initialize_with_mode" is inconsistent with mode_compute=0 or empty mode_file.')
             else
-                options_.mh_posterior_mode_estimation=0;
+                options_.mh_posterior_mode_estimation=false;
             end
         end
 
@@ -335,7 +335,7 @@ if init
         end
 
         % moreover slice must be associated to:
-        %     options_.mh_posterior_mode_estimation = 0;
+        %     options_.mh_posterior_mode_estimation = false;
         % this is done below, but perhaps preprocessing should do this?
 
         if ~isempty(posterior_sampler_options.mode)
@@ -425,7 +425,7 @@ if strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
     end
     % needs to be re-set to zero otherwise posterior analysis is filtered
     % out in dynare_estimation_1.m
-    options_.mh_posterior_mode_estimation = 0;
+    options_.mh_posterior_mode_estimation = false;
 end
 
 return
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 0f91e457e9b00b16b948d206eed948f927379f26..2424473e03cd88794e70688b8e5442565cc4e916 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -402,7 +402,6 @@ options_.mh_initialize_from_previous_mcmc.status = false;
 options_.mh_initialize_from_previous_mcmc.directory = '';
 options_.mh_initialize_from_previous_mcmc.record = '';
 options_.mh_initialize_from_previous_mcmc.prior = '';
-options_.mh_mode = 1;
 options_.mh_nblck = 2;
 options_.mh_recover = false;
 options_.mh_replic = 20000;
@@ -427,7 +426,7 @@ options_.moments_varendo = false;
 options_.nk = 1;
 options_.noconstant = false;
 options_.nodiagnostic = false;
-options_.mh_posterior_mode_estimation = 0;
+options_.mh_posterior_mode_estimation = false;
 options_.prefilter = 0;
 options_.presample = 0;
 options_.prior_trunc = 1e-10;
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 32160ad58ea1ae312bace2bcd7dae00505ae1916..49df96c69e4c3a344c4d37b8025791d34debcdb8 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -171,7 +171,7 @@ catch % if check fails, provide info on using calibration if present
     rethrow(e);
 end
 
-if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0
+if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation
     if options_.order==1 && ~options_.particle.status
         if options_.smoother
             if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter
@@ -475,6 +475,8 @@ if (any(bayestopt_.pshape  >0 ) && options_.mh_replic) || ...
     %% Here I discard first mh_drop percent of the draws:
     CutSample(M_, options_, estim_params_);
     if options_.mh_posterior_mode_estimation
+        [~,~,posterior_mode,~] = compute_mh_covariance_matrix();
+        oo_=fill_mh_mode(posterior_mode',NaN(length(posterior_mode),1),M_,options_,estim_params_,bayestopt_,oo_,'posterior');
         %reset qz_criterium
         options_.qz_criterium=qz_criterium_old;
         return
diff --git a/matlab/fill_mh_mode.m b/matlab/fill_mh_mode.m
new file mode 100644
index 0000000000000000000000000000000000000000..5682e1f7fa2814047b9f6d36d3d417d0d05aa8a1
--- /dev/null
+++ b/matlab/fill_mh_mode.m
@@ -0,0 +1,93 @@
+function oo_=fill_mh_mode(xparam1,stdh,M_,options_,estim_params_,bayestopt_,oo_, field_name)
+%function oo_=fill_mh_mode(xparam1,stdh,M_,options_,estim_params_,bayestopt_,oo_, field_name)
+%
+% INPUTS
+%   o xparam1       [double]   (p*1) vector of estimate parameters.
+%   o stdh          [double]   (p*1) vector of estimate parameters.
+%   o M_                        Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
+%   o estim_params_             Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
+%   o options_                  Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
+%   o bayestopt_                Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}).
+%   o oo_                       Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
+%
+% OUTPUTS
+%   o oo_                       Matlab's structure gathering the results
+%
+% SPECIAL REQUIREMENTS
+%   None.
+
+% Copyright (C) 2005-2021 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 <https://www.gnu.org/licenses/>.
+
+nvx = estim_params_.nvx;  % Variance of the structural innovations (number of parameters).
+nvn = estim_params_.nvn;  % Variance of the measurement innovations (number of parameters).
+ncx = estim_params_.ncx;  % Covariance of the structural innovations (number of parameters).
+ncn = estim_params_.ncn;  % Covariance of the measurement innovations (number of parameters).
+np  = estim_params_.np ;  % Number of deep parameters.
+
+if np
+    ip = nvx+nvn+ncx+ncn+1;
+    for i=1:np
+        name = bayestopt_.name{ip};
+        oo_.([field_name '_mode']).parameters.(name) = xparam1(ip);
+        oo_.([field_name '_std_at_mode']).parameters.(name) = stdh(ip);
+        ip = ip+1;
+    end
+end
+if nvx
+    ip = 1;
+    for i=1:nvx
+        k = estim_params_.var_exo(i,1);
+        name = M_.exo_names{k};
+        oo_.([field_name '_mode']).shocks_std.(name)= xparam1(ip);
+        oo_.([field_name '_std_at_mode']).shocks_std.(name) = stdh(ip);
+        ip = ip+1;
+    end
+end
+if nvn
+    ip = nvx+1;
+    for i=1:nvn
+        name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
+        oo_.([field_name '_mode']).measurement_errors_std.(name) = xparam1(ip);
+        oo_.([field_name '_std_at_mode']).measurement_errors_std.(name) = stdh(ip);
+        ip = ip+1;
+    end
+end
+
+if ncx
+    ip = nvx+nvn+1;
+    for i=1:ncx
+        k1 = estim_params_.corrx(i,1);
+        k2 = estim_params_.corrx(i,2);
+        NAME = [M_.exo_names{k1} '_' M_.exo_names{k2}];
+        oo_.([field_name '_mode']).shocks_corr.(name) = xparam1(ip);
+        oo_.([field_name '_std_at_mode']).shocks_corr.(name) = stdh(ip);
+        ip = ip+1;
+    end
+end
+
+if ncn
+    ip = nvx+nvn+ncx+1;
+    for i=1:ncn
+        k1 = estim_params_.corrn(i,1);
+        k2 = estim_params_.corrn(i,2);
+        NAME = [M_.endo_names{k1} '_' M_.endo_names{k2}];
+        oo_.([field_name '_mode']).measurement_errors_corr.(name) = xparam1(ip);
+        oo_.([field_name '_std_at_mode']).measurement_errors_corr.(name) = stdh(ip);
+        ip = ip+1;
+    end
+end
\ No newline at end of file
diff --git a/matlab/marginal_density.m b/matlab/marginal_density.m
index bd10e86ea0a4df50dd339c097613c3fcf34a811c..6cd7e816b54d59a385ff47b1e3d4f754cfbf46af 100644
--- a/matlab/marginal_density.m
+++ b/matlab/marginal_density.m
@@ -126,85 +126,3 @@ end
 
 oo_.MarginalDensity.ModifiedHarmonicMean = mean(marginal(:,2));
 
-return
-
-function oo_=fill_mh_mode(xparam1,stdh,M_,options_,estim_params_,bayestopt_,oo_, field_name)
-%function oo_=fill_mh_mode(xparam1,stdh,M_,options_,estim_params_,bayestopt_,oo_, field_name)
-%
-% INPUTS
-%   o xparam1       [double]   (p*1) vector of estimate parameters.
-%   o stdh          [double]   (p*1) vector of estimate parameters.
-%   o M_                        Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
-%   o estim_params_             Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
-%   o options_                  Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
-%   o bayestopt_                Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}).
-%   o oo_                       Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
-%
-% OUTPUTS
-%   o oo_                       Matlab's structure gathering the results
-%
-% SPECIAL REQUIREMENTS
-%   None.
-
-
-nvx = estim_params_.nvx;  % Variance of the structural innovations (number of parameters).
-nvn = estim_params_.nvn;  % Variance of the measurement innovations (number of parameters).
-ncx = estim_params_.ncx;  % Covariance of the structural innovations (number of parameters).
-ncn = estim_params_.ncn;  % Covariance of the measurement innovations (number of parameters).
-np  = estim_params_.np ;  % Number of deep parameters.
-nx  = nvx+nvn+ncx+ncn+np; % Total number of parameters to be estimated.
-
-if np
-    ip = nvx+nvn+ncx+ncn+1;
-    for i=1:np
-        name = bayestopt_.name{ip};
-        eval(['oo_.' field_name '_mode.parameters.' name ' = xparam1(ip);']);
-        eval(['oo_.' field_name '_std_at_mode.parameters.' name ' = stdh(ip);']);
-        ip = ip+1;
-    end
-end
-if nvx
-    ip = 1;
-    for i=1:nvx
-        k = estim_params_.var_exo(i,1);
-        name = M_.exo_names{k};
-        eval(['oo_.' field_name '_mode.shocks_std.' name ' = xparam1(ip);']);
-        eval(['oo_.' field_name '_std_at_mode.shocks_std.' name ' = stdh(ip);']);
-        ip = ip+1;
-    end
-end
-if nvn
-    ip = nvx+1;
-    for i=1:nvn
-        name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
-        eval(['oo_.' field_name '_mode.measurement_errors_std.' name ' = xparam1(ip);']);
-        eval(['oo_.' field_name '_std_at_mode.measurement_errors_std.' name ' = stdh(ip);']);
-        ip = ip+1;
-    end
-end
-
-if ncx
-    ip = nvx+nvn+1;
-    for i=1:ncx
-        k1 = estim_params_.corrx(i,1);
-        k2 = estim_params_.corrx(i,2);
-        NAME = [M_.exo_names{k1} '_' M_.exo_names{k2}];
-        eval(['oo_.' field_name '_mode.shocks_corr.' NAME ' = xparam1(ip);']);
-        eval(['oo_.' field_name '_std_at_mode.shocks_corr.' NAME ' = stdh(ip);']);
-        ip = ip+1;
-    end
-end
-
-if ncn
-    ip = nvx+nvn+ncx+1;
-    for i=1:ncn
-        k1 = estim_params_.corrn(i,1);
-        k2 = estim_params_.corrn(i,2);
-        NAME = [M_.endo_names{k1} '_' M_.endo_names{k2}];
-        eval(['oo_.' field_name '_mode.measurement_errors_corr.' NAME ' = xparam1(ip);']);
-        eval(['oo_.' field_name '_std_at_mode.measurement_errors_corr.' NAME ' = stdh(ip);']);
-        ip = ip+1;
-    end
-end
-
-return
\ No newline at end of file