diff --git a/doc/dynare.texi b/doc/dynare.texi
index 0f01b3fa109650866205b591e4560f9a73a7691c..0d30190babfc9710c4b928f014a0f97cf8e27d50 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -4461,6 +4461,9 @@ Median of the posterior distribution
 @item Std
 Standard deviation of the posterior distribution
 
+@item Variance
+Variance of the posterior distribution
+
 @item deciles
 Deciles of the distribution.
 
@@ -4649,6 +4652,21 @@ oo_.posterior_mean.shocks_std.ex
 oo_.posterior_hpdsup.measurement_errors_corr.gdp_conso
 @end example
 
+@defvr {MATLAB/Octave variable} oo_.MC_record.Seeds
+Variable set by the @code{estimation} command. Stores seeds used in MCMC chains
+@end defvr
+
+@defvr {MATLAB/Octave variable} oo_.MC_record.AcceptationRates
+Variable set by the @code{estimation} command. Stores acceptation rates of the MCMC chains
+@end defvr
+
+@defvr {MATLAB/Octave variable} oo_.MC_record.LastParameters
+Variable set by the @code{estimation} command. Stores parameter vector of final MCMC chain draw
+@end defvr
+
+@defvr {MATLAB/Octave variable} oo_.MC_record.LastLogPost
+Variable set by the @code{estimation} command. Stores log-posterior of final MCMC chain draw
+@end defvr
 
 @deffn Command model_comparison @var{FILENAME}[(@var{DOUBLE})]@dots{};
 @deffnx Command model_comparison (marginal_density = laplace | modifiedharmonicmean) @var{FILENAME}[(@var{DOUBLE})]@dots{};
@@ -5316,7 +5334,7 @@ The discussion of the methodologies and their application is described in
 
 With respect to the previous version of the toolbox, in order to work
 properly, the GSA toolbox no longer requires that the Dynare
-estimation environment is setup.
+estimation environment is set up.
 
 Sensitivity analysis results are saved locally in @code{<mod_file>/GSA},
 where @code{<mod_file>.mod} is the name of the DYNARE model file.
diff --git a/matlab/adaptive_metropolis_hastings.m b/matlab/adaptive_metropolis_hastings.m
index ea71612f06cf4f12f51a595a45e687d09a3fd3ef..2b0cb451a83745af7397bb2b8537d87f244f4fd4 100644
--- a/matlab/adaptive_metropolis_hastings.m
+++ b/matlab/adaptive_metropolis_hastings.m
@@ -1,4 +1,4 @@
-function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
+function record=adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
 %function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
 % Random walk Metropolis-Hastings algorithm. 
 % 
@@ -11,7 +11,7 @@ function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds
 %   o varargin              list of argument following mh_bounds
 %  
 % OUTPUTS 
-%   None  
+%   o record     [struct]   structure describing the iterations
 %
 % ALGORITHM 
 %   Metropolis-Hastings.       
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 45e47397f55ea11d718d04b3afda4bdd3584d2f2..f4a6dcadadc9d800fe3f50027425fe0baf9a20fd 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -156,7 +156,7 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.m
         oo_.Smoother.SteadyState = ys;
         oo_.Smoother.TrendCoeffs = trend_coeff;
         if options_.filter_covariance
-            oo_.Smoother.variance = P;
+            oo_.Smoother.Variance = P;
         end
         i_endo = bayestopt_.smoother_saved_var_list;
         if options_.nk ~= 0
@@ -482,12 +482,12 @@ if options_.mode_check == 1 && ~options_.mh_posterior_mode_estimation
 end
 
 oo_.posterior.optimization.mode = xparam1;
-oo_.posterior.optimization.variance = [];
+oo_.posterior.optimization.Variance = [];
 if ~options_.mh_posterior_mode_estimation
     if options_.cova_compute
         invhess = inv(hh);
         stdh = sqrt(diag(invhess));
-        oo_.posterior.optimization.variance = invhess;
+        oo_.posterior.optimization.Variance = invhess;
     end
 else
     variances = bayestopt_.p2.*bayestopt_.p2;
@@ -921,7 +921,7 @@ if (any(bayestopt_.pshape  >0 ) && options_.mh_replic) || ...
         ana_deriv = options_.analytic_derivation;
         options_.analytic_derivation = 0;
         if options_.cova_compute
-            feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+            oo_.MC_record=feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
         else
             error('I Cannot start the MCMC because the hessian of the posterior kernel at the mode was not computed.')
         end
@@ -943,7 +943,7 @@ if (any(bayestopt_.pshape  >0 ) && options_.mh_replic) || ...
             if ~options_.nograph
                 oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt_, oo_);
             end
-            [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.variance] ...
+            [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.Variance] ...
                 = GetPosteriorMeanVariance(M_,options_.mh_drop);
         else
             load([M_.fname '_results'],'oo_');
@@ -974,7 +974,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
     [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,dataset_.info.ntobs,dataset_.data,dataset_.missing.aindex,dataset_.missing.state);
     oo_.Smoother.SteadyState = ys;
     oo_.Smoother.TrendCoeffs = trend_coeff;
-    oo_.Smoother.variance = P;
+    oo_.Smoother.Variance = P;
     i_endo = bayestopt_.smoother_saved_var_list;
     if options_.nk ~= 0
         oo_.FilteredVariablesKStepAhead = aK(options_.filter_step_ahead, ...
diff --git a/matlab/independent_metropolis_hastings.m b/matlab/independent_metropolis_hastings.m
index bf198383fc898f8f697bcbede5f908ee829eb7e3..8e24ad1788fc66afb17bf867a90ebbf655e16e17 100644
--- a/matlab/independent_metropolis_hastings.m
+++ b/matlab/independent_metropolis_hastings.m
@@ -1,4 +1,4 @@
-function independent_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
+function record=independent_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
 
 % Independent Metropolis-Hastings algorithm. 
 % 
@@ -11,7 +11,7 @@ function independent_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bou
 %   o varargin              list of argument following mh_bounds
 %  
 % OUTPUTS 
-%   None  
+%   o record     [struct]   structure describing the iterations
 %
 % ALGORITHM 
 %   Metropolis-Hastings.       
@@ -103,7 +103,7 @@ else
     [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'independent_metropolis_hastings_core', localVars, globalVars, options_.parallel_info);
     for j=1:totCPU,
         offset = sum(nBlockPerCPU(1:j-1))+fblck-1;
-        record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)));
+        record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)));
         record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:)=fout(j).record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:);
         record.AcceptationRates(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.AcceptationRates(offset+1:sum(nBlockPerCPU(1:j)));
         record.Seeds(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.Seeds(offset+1:sum(nBlockPerCPU(1:j)));
diff --git a/matlab/independent_metropolis_hastings_core.m b/matlab/independent_metropolis_hastings_core.m
index 04b2683aa0a76c75122ca128a161eb5b06652a4a..efb50965e05ee5ad386b23afc057e940617f9409 100644
--- a/matlab/independent_metropolis_hastings_core.m
+++ b/matlab/independent_metropolis_hastings_core.m
@@ -229,7 +229,7 @@ for b = fblck:nblck,
             jsux = 0;
             if j == nruns(b) % I record the last draw...
                 record.LastParameters(b,:) = x2(end,:);
-                record.LastLogLiK(b) = logpo2(end);
+                record.LastLogPost(b) = logpo2(end);
             end
             % size of next file in chain b
             InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);
diff --git a/matlab/metropolis_hastings_initialization.m b/matlab/metropolis_hastings_initialization.m
index 06dc8a4ca9fcbadfb16a9a9153fa031c5abc2f55..e56ddcf5f40ef2827d59067f6afd688f3d58546e 100644
--- a/matlab/metropolis_hastings_initialization.m
+++ b/matlab/metropolis_hastings_initialization.m
@@ -194,7 +194,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
     record.InitialParameters = ix2;
     record.InitialLogLiK = ilogpo2;
     record.LastParameters = zeros(nblck,npar);
-    record.LastLogLiK = zeros(nblck,1);
+    record.LastLogPost = zeros(nblck,1);
     record.LastFileNumber = AnticipatedNumberOfFiles ;
     record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
     save([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');  
@@ -253,7 +253,7 @@ elseif options_.load_mh_file && ~options_.mh_recover
     else
         NewFile = ones(nblck,1)*(LastFileNumber+1);
     end
-    ilogpo2 = record.LastLogLiK;
+    ilogpo2 = record.LastLogPost;
     ix2 = record.LastParameters;
     fblck = 1;
     fline = ones(nblck,1)*(LastLineNumber+1);
@@ -292,7 +292,7 @@ elseif options_.mh_recover
     end
     %% Default initialization:
     if OldMh
-        ilogpo2 = record.LastLogLiK;
+        ilogpo2 = record.LastLogPost;
         ix2 = record.LastParameters;
     else
         ilogpo2 = record.InitialLogLiK;
diff --git a/matlab/random_walk_metropolis_hastings.m b/matlab/random_walk_metropolis_hastings.m
index 453e4d615939ad960bddde01aa23d17d4372afb2..1c75e210290b385b1dd7195c74c00721efb5060b 100644
--- a/matlab/random_walk_metropolis_hastings.m
+++ b/matlab/random_walk_metropolis_hastings.m
@@ -148,7 +148,7 @@ else
     [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'random_walk_metropolis_hastings_core', localVars, globalVars, options_.parallel_info);
     for j=1:totCPU,
         offset = sum(nBlockPerCPU(1:j-1))+fblck-1;
-        record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)));
+        record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)));
         record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:)=fout(j).record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:);
         record.AcceptationRates(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.AcceptationRates(offset+1:sum(nBlockPerCPU(1:j)));
         record.Seeds(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.Seeds(offset+1:sum(nBlockPerCPU(1:j)));
diff --git a/matlab/random_walk_metropolis_hastings_core.m b/matlab/random_walk_metropolis_hastings_core.m
index 69062f66c735746eb316f0d69b9b9f360e057342..a92c3c21c72f1071f04197f4d3c72033f854bdef 100644
--- a/matlab/random_walk_metropolis_hastings_core.m
+++ b/matlab/random_walk_metropolis_hastings_core.m
@@ -247,7 +247,7 @@ for b = fblck:nblck,
             jsux = 0;
             if j == nruns(b) % I record the last draw...
                 record.LastParameters(b,:) = x2(end,:);
-                record.LastLogLiK(b) = logpo2(end);
+                record.LastLogPost(b) = logpo2(end);
             end
             % size of next file in chain b
             InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);