Merge pull request #1297 from JohannesPfeifer/estimation_fixes

Bugfixes and cosmetic changes related to estimation

matlab/GetPosteriorMeanVariance.m

 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

matlab/dyn_autocorr.m

-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.
 %

matlab/posterior_sampler_core.m

@@ -126,7 +126,11 @@ 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)
+        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);            
     catch

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