diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m
index 4932e304262a4594c202a30a65112ecc55240014..01db632a8afb62ad4e798f2c97de0ac47bb58e15 100644
--- a/matlab/dynare_estimation_init.m
+++ b/matlab/dynare_estimation_init.m
@@ -500,18 +500,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