diff --git a/matlab/GetPosteriorMeanVariance.m b/matlab/GetPosteriorMeanVariance.m
index e364180917773afc7730d116b748aeca6c7047f1..d1e2935873b0e5a969446e4c2572c31756469140 100644
--- a/matlab/GetPosteriorMeanVariance.m
+++ b/matlab/GetPosteriorMeanVariance.m
@@ -1,6 +1,6 @@
 function [mean,variance] = GetPosteriorMeanVariance(M,drop)
 
-% Copyright (C) 2012, 2013 Dynare Team
+% Copyright (C) 2012-2016 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -26,29 +26,28 @@ function [mean,variance] = GetPosteriorMeanVariance(M,drop)
     NbrBlocks = record.Nblck;
     mean = 0;
     variance = 0;
-    z = [];
     
-    nkept = 0;
+    NbrKeptDraws = 0;
     for i=1:NbrBlocks
-        n = 0;
+        NbrDrawsCurrentBlock = 0;
         for j=1:NbrFiles
-            o = load([BaseName '_mh' int2str(j) '_blck' int2str(i)]);
-            m = size(o.x2,1);
-            if n + m < drop*NbrDraws
-                n = n + m;
+            o = load([BaseName '_mh' int2str(j) '_blck' int2str(i),'.mat']);
+            NbrDrawsCurrentFile = size(o.x2,1);
+            if NbrDrawsCurrentBlock + NbrDrawsCurrentFile <= drop*NbrDraws
+                NbrDrawsCurrentBlock = NbrDrawsCurrentBlock + NbrDrawsCurrentFile;
                 continue
-            elseif n < drop*NbrDraws
-                k = ceil(drop*NbrDraws - n + 1);
-                x2 = o.x2(k:end,:);
+            elseif NbrDrawsCurrentBlock < drop*NbrDraws
+                FirstDraw = ceil(drop*NbrDraws - NbrDrawsCurrentBlock + 1);
+                x2 = o.x2(FirstDraw:end,:);
             else
                 x2 = o.x2;
             end
-            z =[z; x2];        
-            p = size(x2,1);
-            mean = (nkept*mean + sum(x2)')/(nkept+p);
-            x = bsxfun(@minus,x2,mean');
-            variance = (nkept*variance + x'*x)/(nkept+p);
-            n = n + m;
-            nkept = nkept + p;
+            NbrKeptDrawsCurrentFile = size(x2,1);
+            %recursively compute mean and variance
+            mean = (NbrKeptDraws*mean + sum(x2)')/(NbrKeptDraws+NbrKeptDrawsCurrentFile);
+            x2Demeaned = bsxfun(@minus,x2,mean');
+            variance = (NbrKeptDraws*variance + x2Demeaned'*x2Demeaned)/(NbrKeptDraws+NbrKeptDrawsCurrentFile);
+            NbrDrawsCurrentBlock = NbrDrawsCurrentBlock + NbrDrawsCurrentFile;
+            NbrKeptDraws = NbrKeptDraws + NbrKeptDrawsCurrentFile;
         end
     end
diff --git a/matlab/dyn_autocorr.m b/matlab/dyn_autocorr.m
index a975511334acae02a1f5666e85fe2273fd5f9d3d..43534d7bee736fb575aa6b60c76e7505de6968fc 100644
--- a/matlab/dyn_autocorr.m
+++ b/matlab/dyn_autocorr.m
@@ -1,5 +1,5 @@
-function acf = autocorr(y, ar)
-% function acf = autocorr(y, ar)
+function acf = dyn_autocorr(y, ar)
+% function acf = dyn_autocorr(y, ar)
 % autocorrelation function of y
 %
 % INPUTS
@@ -12,7 +12,7 @@ function acf = autocorr(y, ar)
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright (C) 2015 Dynare Team
+% Copyright (C) 2015-16 Dynare Team
 %
 % This file is part of Dynare.
 %
diff --git a/matlab/posterior_sampler_core.m b/matlab/posterior_sampler_core.m
index c677e02cceb28627103039868b0ee1948ff569e0..d88ae742ac633dd75ae2c12bbadc7169a2d93f08 100644
--- a/matlab/posterior_sampler_core.m
+++ b/matlab/posterior_sampler_core.m
@@ -126,9 +126,13 @@ for curr_block = fblck:nblck,
         % Matlab/Octave cluster). Therefore the trap.
         %
         % Set the random number generator type (the seed is useless but needed by the function)
-        set_dynare_seed(options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed);
+        if ~isoctave
+            set_dynare_seed(options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed);
+        else
+            set_dynare_seed(options_.DynareRandomStreams.seed+curr_block);
+        end
         % Set the state of the RNG
-        set_dynare_random_generator_state(record.InitialSeeds(curr_block).Unifor, record.InitialSeeds(curr_block).Normal);
+        set_dynare_random_generator_state(record.InitialSeeds(curr_block).Unifor, record.InitialSeeds(curr_block).Normal);            
     catch
         % If the state set by master is incompatible with the slave, we only reseed 
         set_dynare_seed(options_.DynareRandomStreams.seed+curr_block);
diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m
index 3adca10a50105e8ae45dfdb6c0ab13bfe47b19c1..07af053851a3a10dc6fb6dce6c1186e93589e88c 100644
--- a/matlab/posterior_sampler_initialization.m
+++ b/matlab/posterior_sampler_initialization.m
@@ -458,7 +458,7 @@ function [d,bayestopt_]=set_proposal_density_to_previous_value(record,options_,b
                 error(fprintf('Estimation::mcmc: The posterior_sampling_method differs from the one of the original chain. Please reset it to %s',record.MCMC_sampler))
             end 
         end
-        fprintf('Estimation::mcmc: Recovering the previous proposal density\n.')
+        fprintf('Estimation::mcmc: Recovering the previous proposal density.\n')
         d=record.ProposalCovariance;
         bayestopt_.jscale=record.ProposalScaleVec;
     else