diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m
index ae123f832..a6696e944 100644
--- a/matlab/disp_th_moments.m
+++ b/matlab/disp_th_moments.m
@@ -38,7 +38,8 @@ function disp_th_moments(dr,var_list)
   [oo_.gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_,options_);
   m = dr.ys(ivar);
   non_stationary_vars = setdiff(1:length(ivar),stationary_vars);
-  m(ivar(non_stationary_vars)) = NaN;
+  ivar1 = intersect(non_stationary_vars,ivar);
+  m(ivar1) = NaN;
 
   
   i1 = find(abs(diag(oo_.gamma_y{1})) > 1e-12);
diff --git a/matlab/get_variance_of_endogenous_variables.m b/matlab/get_variance_of_endogenous_variables.m
index 21542d2c4..77db2c5f7 100644
--- a/matlab/get_variance_of_endogenous_variables.m
+++ b/matlab/get_variance_of_endogenous_variables.m
@@ -33,6 +33,8 @@ function [vx1,i_ns] = get_variance_of_endogenous_variables(dr,i_var)
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
   global M_ options_
+   
+  endo_nbr = M_.endo_nbr;
   
   Sigma_e = M_.Sigma_e;
   
@@ -48,6 +50,7 @@ function [vx1,i_ns] = get_variance_of_endogenous_variables(dr,i_var)
   [vx,u] = lyapunov_symm(A,B*Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
   
   if size(u,2) > 0
+    i_stat_0 = find(any(abs(A*u) < options_.Schur_vec_tol,2));
     i_stat = find(any(abs(ghx*u) < options_.Schur_vec_tol,2));
   
     ghx = ghx(i_stat,:);
@@ -57,5 +60,5 @@ function [vx1,i_ns] = get_variance_of_endogenous_variables(dr,i_var)
   end
   
   vx1 = Inf*ones(n,n);
-  vx1(i_stat,i_stat) = ghx*vx*ghx'+ghu*Sigma_e*ghu';
-  
\ No newline at end of file
+  vx1(i_stat,i_stat) = ghx(:,i_stat_0)*vx(i_stat_0,i_stat_0)*ghx(:,i_stat_0)'+ghu*Sigma_e*ghu';
+