diff --git a/matlab/distributions/inverse_gamma_specification.m b/matlab/distributions/inverse_gamma_specification.m
index e5136381ad7591d55a029d6a25b924eb6b08a74f..e5ca2d3ba2134fb202dc246cc2a146d91243c5ca 100644
--- a/matlab/distributions/inverse_gamma_specification.m
+++ b/matlab/distributions/inverse_gamma_specification.m
@@ -85,7 +85,7 @@ s = [];
 nu = [];
 
 sigma = sqrt(sigma2);
-mu2 = mu*mu;
+mu2 = (mu-lb)*(mu-lb);
 
 if type == 2       % Inverse Gamma 2
     nu   = 2*(2+mu2/sigma2);
@@ -132,10 +132,10 @@ elseif type == 1   % Inverse Gamma 1
         end
         s = (sigma2+mu2)*(nu-2);
         if check_solution_flag
-            if abs(log(mu)-log(sqrt(s/2))-gammaln((nu-1)/2)+gammaln(nu/2))>1e-7
+            if abs(log(mu-lb)-log(sqrt(s/2))-gammaln((nu-1)/2)+gammaln(nu/2))>1e-7
                 error('inverse_gamma_specification:: Failed in solving for the hyperparameters!');
             end
-            if abs(sigma-sqrt(s/(nu-2)-mu*mu))>1e-7
+            if abs(sigma-sqrt(s/(nu-2)-(mu-lb)*(mu-lb)))>1e-7
                 error('inverse_gamma_specification:: Failed in solving for the hyperparameters!');
             end
         end
diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m
index 76fa750ffb741bfc79a6cb5e290409889e70719a..8238f7a9a4e2cedc713f0b49fc167ce0f67f622f 100644
--- a/matlab/dynare_estimation_init.m
+++ b/matlab/dynare_estimation_init.m
@@ -498,18 +498,18 @@ if options_.heteroskedastic_filter
             'for heteroskedastic_filter'])
     end
 
-    M_.heteroskedastic_shocks.Qvalue = NaN(M_.exo_nbr,options_.nobs);
-    M_.heteroskedastic_shocks.Qscale = NaN(M_.exo_nbr,options_.nobs);
+    M_.heteroskedastic_shocks.Qvalue = NaN(M_.exo_nbr,options_.nobs+1);
+    M_.heteroskedastic_shocks.Qscale = NaN(M_.exo_nbr,options_.nobs+1);
 
     for k=1:length(M_.heteroskedastic_shocks.Qvalue_orig)
         v = M_.heteroskedastic_shocks.Qvalue_orig(k);
-        temp_periods=v.periods(v.periods<options_.nobs+options_.first_obs);
+        temp_periods=v.periods(v.periods<options_.nobs+options_.first_obs+1);
         temp_periods=temp_periods(temp_periods>=options_.first_obs);
         M_.heteroskedastic_shocks.Qvalue(v.exo_id, temp_periods-(options_.first_obs-1)) = v.value^2;
     end
     for k=1:length(M_.heteroskedastic_shocks.Qscale_orig)
         v = M_.heteroskedastic_shocks.Qscale_orig(k);
-        temp_periods=v.periods(v.periods<options_.nobs+options_.first_obs);
+        temp_periods=v.periods(v.periods<options_.nobs+options_.first_obs+1);
         temp_periods=temp_periods(temp_periods>=options_.first_obs);
         M_.heteroskedastic_shocks.Qscale(v.exo_id, temp_periods-(options_.first_obs-1)) = v.scale^2;
     end
diff --git a/matlab/kalman/get_Qvec_heteroskedastic_filter.m b/matlab/kalman/get_Qvec_heteroskedastic_filter.m
index 7946b48d51b5744c14873846626992fbe70a5fe6..ae4e6886895a2b35dfbe22716dec563ba0030fb4 100644
--- a/matlab/kalman/get_Qvec_heteroskedastic_filter.m
+++ b/matlab/kalman/get_Qvec_heteroskedastic_filter.m
@@ -27,7 +27,7 @@ function Qvec=get_Qvec_heteroskedastic_filter(Q,smpl,M_)
 
 isqdiag = all(all(abs(Q-diag(diag(Q)))<1e-14)); % ie, the covariance matrix is diagonal...
 Qvec=repmat(Q,[1 1 smpl+1]);
-for k=1:smpl
+for k=1:smpl+1
     inx = ~isnan(M_.heteroskedastic_shocks.Qvalue(:,k));
     if any(inx)
         if isqdiag
diff --git a/matlab/posterior_sampler_core.m b/matlab/posterior_sampler_core.m
index 551ba4dafc1be9e92928fca572fa647ea8ed337e..fb4bf8f992a4212bbd772784770015b19f906faf 100644
--- a/matlab/posterior_sampler_core.m
+++ b/matlab/posterior_sampler_core.m
@@ -141,42 +141,51 @@ for curr_block = fblck:nblck
         options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.seed+curr_block);
     end
     mh_recover_flag=0;
-    if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood
-        load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'])
-        x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)];
-        logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)];
+    if options_.mh_recover && exist([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'file')==2 && OpenOldFile(curr_block)
+        % this should be done whatever value of load_mh_file
+        load([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']);
+        draw_iter = size(neval_this_chain,2)+1;
+        draw_index_current_file = draw_iter+fline(curr_block)-1;
+        feval_this_chain = sum(sum(neval_this_chain));
+        feval_this_file = sum(sum(neval_this_chain));
+        if feval_this_chain>draw_index_current_file-fline(curr_block)
+            % non Metropolis type of sampler
+            accepted_draws_this_chain = draw_index_current_file-fline(curr_block);
+            accepted_draws_this_file = draw_index_current_file-fline(curr_block);
+        else
+            accepted_draws_this_chain = 0;
+            accepted_draws_this_file = 0;
+        end
+        mh_recover_flag=1;
+        set_dynare_random_generator_state(LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal);
+        last_draw(curr_block,:)=x2(draw_index_current_file-1,:);
+        last_posterior(curr_block)=logpo2(draw_index_current_file-1);
         OpenOldFile(curr_block) = 0;
     else
-        if options_.mh_recover && exist([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'file')==2
-            load([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']);
-            draw_iter = size(neval_this_chain,2)+1;
-            draw_index_current_file = draw_iter;
-            feval_this_chain = sum(sum(neval_this_chain));
-            feval_this_file = sum(sum(neval_this_chain));
-            if feval_this_chain>draw_iter-1
-                % non Metropolis type of sampler
-                accepted_draws_this_chain = draw_iter-1;
-                accepted_draws_this_file = draw_iter-1;
-            else
-                accepted_draws_this_chain = 0;
-                accepted_draws_this_file = 0;
-            end
-            mh_recover_flag=1;
-            set_dynare_random_generator_state(LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal);
-            last_draw(curr_block,:)=x2(draw_iter-1,:);
-            last_posterior(curr_block)=logpo2(draw_iter-1);
-
+        if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood
+            load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'])
+            x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)];
+            logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)];
+            OpenOldFile(curr_block) = 0;
         else
 
             x2 = zeros(InitSizeArray(curr_block),npar);
             logpo2 = zeros(InitSizeArray(curr_block),1);
         end
     end
+    if mh_recover_flag==0
+        accepted_draws_this_chain = 0;
+        accepted_draws_this_file = 0;
+        feval_this_chain = 0;
+        feval_this_file = 0;
+        draw_iter = 1;
+        draw_index_current_file = fline(curr_block); %get location of first draw in current block
+    end
     %Prepare waiting bars
     if whoiam
         refresh_rate = sampler_options.parallel_bar_refresh_rate;
         bar_title = sampler_options.parallel_bar_title;
-        prc0=(curr_block-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode);
+        prc0=(curr_block-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode)+(draw_iter-1)/nruns(curr_block);
         hh_fig = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
     else
         refresh_rate = sampler_options.serial_bar_refresh_rate;
@@ -184,14 +193,6 @@ for curr_block = fblck:nblck
         hh_fig = dyn_waitbar(0,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
         set(hh_fig,'Name',bar_title);
     end
-    if mh_recover_flag==0
-        accepted_draws_this_chain = 0;
-        accepted_draws_this_file = 0;
-        feval_this_chain = 0;
-        feval_this_file = 0;
-        draw_iter = 1;
-        draw_index_current_file = fline(curr_block); %get location of first draw in current block
-    end
     sampler_options.curr_block = curr_block;
     while draw_iter <= nruns(curr_block)
 
diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m
index 595715a2bc53dae2010d8ab4959d23f65724a577..3ed18037f1e0b34324cbccb4149450eb02cfc1c5 100644
--- a/matlab/posterior_sampler_initialization.m
+++ b/matlab/posterior_sampler_initialization.m
@@ -170,6 +170,12 @@ if ~options_.load_mh_file && ~options_.mh_recover
         ilogpo2 = zeros(NumberOfBlocks,1);
         ilogpo2(:,1) = record0.LastLogPost;
         ix2(:,IA) = record0.LastParameters(:,IB);
+        for j=1:NumberOfBlocks
+            if not(all(ix2(j,:)' >= mh_bounds.lb) && all(ix2(j,:)' <= mh_bounds.ub))
+                new_estimated_parameters = logical(new_estimated_parameters + (ix2(j,:)' < mh_bounds.lb));
+                new_estimated_parameters = logical(new_estimated_parameters + (ix2(j,:)' > mh_bounds.ub));
+            end
+        end
     else
         new_estimated_parameters = true(1,npar);
     end
@@ -458,7 +464,7 @@ elseif options_.mh_recover
     AllMhFiles = dir([BaseName '_mh*_blck*.mat']);
     TotalNumberOfMhFiles = length(AllMhFiles)-length(dir([BaseName '_mh_tmp*_blck*.mat']));
     % Quit if no crashed mcmc chain can be found as there are as many files as expected
-    if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles)
+    if (ExpectedNumberOfMhFilesPerBlock>LastFileNumberInThePreviousMh) && (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles)
         if isnumeric(options_.parallel)
             fprintf('%s: It appears that you don''t need to use the mh_recover option!\n',dispString);
             fprintf('                  You have to edit the mod file and remove the mh_recover option\n');
@@ -469,15 +475,23 @@ elseif options_.mh_recover
     % 2. Something needs to be done; find out what
     % Count the number of saved mh files per block.
     NumberOfMhFilesPerBlock = zeros(NumberOfBlocks,1);
+    is_chain_complete = true(NumberOfBlocks,1);
     for b = 1:NumberOfBlocks
         BlckMhFiles = dir([BaseName '_mh*_blck' int2str(b) '.mat']);
         NumberOfMhFilesPerBlock(b) = length(BlckMhFiles)-length(dir([BaseName '_mh_tmp*_blck' int2str(b) '.mat']));
+        if (ExpectedNumberOfMhFilesPerBlock==LastFileNumberInThePreviousMh)
+            % here we need to check for the number of lines
+            tmpdata = load([BaseName '_mh' int2str(LastFileNumberInThePreviousMh) '_blck' int2str(b) '.mat'],'logpo2');
+            if length(tmpdata.logpo2) == LastLineNumberInThePreviousMh
+                is_chain_complete(b) = false;
+            end
+        end
     end
     % Find FirstBlock (First block), an integer targeting the crashed mcmc chain.
     FirstBlock = 1; %initialize
     FBlock = zeros(NumberOfBlocks,1);
     while FirstBlock <= NumberOfBlocks
-        if  NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock
+        if  (NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock) || not(is_chain_complete(FirstBlock))
             fprintf('%s: Chain %u is not complete!\n', dispString,FirstBlock);
             FBlock(FirstBlock)=1;
         else