diff --git a/matlab/lyapunov_symm.m b/matlab/lyapunov_symm.m
index 72ae5b086..1c1be1e16 100644
--- a/matlab/lyapunov_symm.m
+++ b/matlab/lyapunov_symm.m
@@ -1,47 +1,61 @@
-% solves x-a*x*a'=b for b (and then x) symmetrical
-function x=lyapunov_symm(a,b)
-  n = size(b,1);
-  if n == 1
-    x=b/(1-a*a);
-    return
-  end
-  x=zeros(n,n);
-  [u,t]=schur(a);
-  b=u'*b*u;
-  for i=n:-1:2
-    if t(i,i-1) == 0
-      if i == n
-	c = zeros(n,1);
-      else
-	c = t(1:i,:)*(x(:,i+1:end)*t(i,i+1:end)')+...
-	    t(i,i)*t(1:i,i+1:end)*x(i+1:end,i);
-      end
-      x(1:i,i)=(eye(i)-t(1:i,1:i)*t(i,i))\(b(1:i,i)+c);
-      x(i,1:i-1)=x(1:i-1,i)';
-    else
-      if i == n
-	c = zeros(n,1);
-	c1 = zeros(n,1);
-      else
-	c = t(1:i,:)*(x(:,i+1:end)*t(i,i+1:end)')+...
-	    t(i,i)*t(1:i,i+1:end)*x(i+1:end,i)+...
-	    t(i,i-1)*t(1:i,i+1:end)*x(i+1:end,i-1);
-	c1 = t(1:i,:)*(x(:,i+1:end)*t(i-1,i+1:end)')+...
-	     t(i-1,i-1)*t(1:i,i+1:end)*x(i+1:end,i-1)+...
-	     t(i-1,i)*t(1:i,i+1:end)*x(i+1:end,i);
-      end
-      z = [eye(i)-t(1:i,1:i)*t(i,i) -t(1:i,1:i)*t(i,i-1);...
-	   -t(1:i,1:i)*t(i-1,i) eye(i)-t(1:i,1:i)*t(i-1,i-1)]...
-	  \[b(1:i,i)+c;b(1:i,i-1)+c1];
-      x(1:i,i) = z(1:i);
-      x(1:i,i-1) = z(i+1:end);
-      x(i,1:i-1)=x(1:i-1,i)';
-      x(i-1,1:i-2)=x(1:i-2,i-1)';
-      i = i - 1;
-    end
-  end
-  if i == 2
-    c = t(1,:)*(x(:,2:end)*t(1,2:end)')+t(1,1)*t(1,2:end)*x(2:end,1);
-    x(1,1)=(b(1,1)+c)/(1-t(1,1)*t(1,1));
-  end
-  x=u*x*u';
\ No newline at end of file
+% solves x-a*x*a'=b for b (and then x) symmetrical
+function [x,ns_var]=lyapunov_symm(a,b)
+  global options_ 
+  
+  info = 0;
+  if size(a,1) == 1
+    x=b/(1-a*a);
+    return
+  end
+  [u,t] = schur(a);
+  if exist('ordeig','builtin')
+    e1 = abs(ordeig(t)) > 2-options_.qz_criterium;
+  else
+    e1 = abs(my_ordeig(t)) > 2-options_.qz_criterium;
+  end
+  k = sum(e1);
+  [u,t] = ordschur(u,t,e1); 
+  n = length(e1)-k;
+  b=u(:,k+1:end)'*b*u(:,k+1:end);
+  t = t(k+1:end,k+1:end);
+  x=zeros(n,n);
+  for i=n:-1:2
+    if t(i,i-1) == 0
+      if i == n
+	c = zeros(n,1);
+      else
+	c = t(1:i,:)*(x(:,i+1:end)*t(i,i+1:end)')+...
+	    t(i,i)*t(1:i,i+1:end)*x(i+1:end,i);
+      end
+      q = eye(i)-t(1:i,1:i)*t(i,i);
+      x(1:i,i) = q\(b(1:i,i)+c);
+      x(i,1:i-1) = x(1:i-1,i)';
+    else
+      if i == n
+	c = zeros(n,1);
+	c1 = zeros(n,1);
+      else
+	c = t(1:i,:)*(x(:,i+1:end)*t(i,i+1:end)')+...
+	    t(i,i)*t(1:i,i+1:end)*x(i+1:end,i)+...
+	    t(i,i-1)*t(1:i,i+1:end)*x(i+1:end,i-1);
+	c1 = t(1:i,:)*(x(:,i+1:end)*t(i-1,i+1:end)')+...
+	     t(i-1,i-1)*t(1:i,i+1:end)*x(i+1:end,i-1)+...
+	     t(i-1,i)*t(1:i,i+1:end)*x(i+1:end,i);
+      end
+      q = [eye(i)-t(1:i,1:i)*t(i,i) -t(1:i,1:i)*t(i,i-1);...
+	   -t(1:i,1:i)*t(i-1,i) eye(i)-t(1:i,1:i)*t(i-1,i-1)];
+      z =  q\[b(1:i,i)+c;b(1:i,i-1)+c1];
+      x(1:i,i) = z(1:i);
+      x(1:i,i-1) = z(i+1:end);
+      x(i,1:i-1)=x(1:i-1,i)';
+      x(i-1,1:i-2)=x(1:i-2,i-1)';
+      i = i - 1;
+    end
+  end
+  if i == 2
+    c = t(1,:)*(x(:,2:end)*t(1,2:end)')+t(1,1)*t(1,2:end)*x(2:end,1);
+    x(1,1)=(b(1,1)+c)/(1-t(1,1)*t(1,1));
+  end
+  x=u(:,k+1:end)*x*u(:,k+1:end)';
+  ns_var = [];
+  ns_var = find(any(abs(u(:,1:k)) > 1e-8,2)); 
diff --git a/matlab/my_ordeig.m b/matlab/my_ordeig.m
new file mode 100644
index 000000000..397f2190d
--- /dev/null
+++ b/matlab/my_ordeig.m
@@ -0,0 +1,18 @@
+function eval = my_ordeig(t)
+  
+  n = size(t,2);
+  eval = zeros(n,1);
+  for i=1:n-1
+    if t(i+1,i) == 0
+      eval(i) = t(i,i);
+    else
+      k = i:i+1;
+      eval(k) = eig(t(k,k));
+      i = i+1;
+    end
+  end
+  if i < n
+    t(n) = t(n,n);
+  end
+  
+      
\ No newline at end of file
diff --git a/matlab/th_autocovariances.m b/matlab/th_autocovariances.m
index a34b746cc..a4d2a6df7 100644
--- a/matlab/th_autocovariances.m
+++ b/matlab/th_autocovariances.m
@@ -1,196 +1,204 @@
-% Copyright (C) 2001 Michel Juillard
-%
-% computes the theoretical auto-covariances, gamma_y, for an AR(p) process 
-% with coefficients dr.ghx and dr.ghu and shock variances M_.Sigma_e
-% for a subset of variables ivar (indices in M_.endo_names)
-% Theoretical HP filtering is available as an option
-
-function gamma_y=th_autocovariances(dr,ivar)
-  global M_ options_
-
-  if sscanf(version('-release'),'%d') < 13
-    warning off
-  else
-    eval('warning off MATLAB:dividebyzero')
-  end
-  nar = options_.ar;
-  gamma_y = cell(nar+1,1);
-  nvar = size(ivar,1);
-  
-  ghx = dr.ghx;
-  ghu = dr.ghu;
-  npred = dr.npred;
-  nstatic = dr.nstatic;
-  kstate = dr.kstate;
-  order = dr.order_var;
-  iv(order) = [1:length(order)];
-  nx = size(ghx,2);
-  
-  ikx = [nstatic+1:nstatic+npred];
-  
-  A = zeros(nx,nx);
-  k0 = kstate(find(kstate(:,2) <= M_.maximum_lag+1),:);
-  i0 = find(k0(:,2) == M_.maximum_lag+1);
-  n0 = length(i0);
-  A(i0,:) = ghx(ikx,:);
-  AS = ghx(:,i0);
-  ghu1 = zeros(nx,M_.exo_nbr);
-  ghu1(i0,:) = ghu(ikx,:);
-  for i=M_.maximum_lag:-1:2
-    i1 = find(k0(:,2) == i);
-    n1 = size(i1,1);
-    j1 = zeros(n1,1);
-    j2 = j1;
-    for k1 = 1:n1
-      j1(k1) = find(k0(1:n0,1)==k0(i1(k1),1));
-      j2(k1) = find(k0(i0,1)==k0(i1(k1),1));
-    end
-    A(i1,i0(j2))=eye(n1);
-    AS(:,j1) = AS(:,j1)+ghx(:,i1);
-    i0 = i1;
-  end
-  b = ghu1*M_.Sigma_e*ghu1';
-  
-  % order of variables with preset variances in ghx and ghu
-  iky = iv(ivar);
-  
-  aa = ghx(iky,:);
-  bb = ghu(iky,:);
-
-  if options_.order == 2
-    %save temp A AS b dr M_.Sigma_e ikx iky nx n0
-    vx =  lyapunov_symm(A,b); 
-    Ex = (dr.ghs2(ikx)+dr.ghxx(ikx,:)*vx(:)+dr.ghuu(ikx,:)*M_.Sigma_e(:))/2;
-    Ex = (eye(n0)-AS(ikx,:))\Ex;
-    gamma_y{nar+3} = AS(iky,:)*Ex+(dr.ghs2(iky)+dr.ghxx(iky,:)*vx(:)+dr.ghuu(iky,:)*M_.Sigma_e(:))/2;
-  end
-  if options_.hp_filter == 0
-    if options_.order < 2
-      vx =  lyapunov_symm(A,b);
-    end
-    gamma_y{1} = aa*vx*aa'+ bb*M_.Sigma_e*bb';
-    k = find(abs(gamma_y{1}) < 1e-12);
-    gamma_y{1}(k) = 0;
-    
-    % autocorrelations
-    if nar > 0
-      vxy = (A*vx*aa'+ghu1*M_.Sigma_e*bb');
-
-      sy = sqrt(diag(gamma_y{1}));
-      sy = sy *sy';
-      gamma_y{2} = aa*vxy./sy;
-      
-      for i=2:nar
-	vxy = A*vxy;
-	gamma_y{i+1} = aa*vxy./sy;
-      end
-    end
-    
-    % variance decomposition
-    if M_.exo_nbr > 1
-      gamma_y{nar+2} = zeros(nvar,M_.exo_nbr);
-      SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord)=M_.Sigma_e+1e-14*eye(M_.exo_nbr);
-      cs = chol(SS)';
-      b1(:,M_.exo_names_orig_ord) = ghu1;
-      b1 = b1*cs;
-      b2(:,M_.exo_names_orig_ord) = ghu(iky,:);
-      b2 = b2*cs;
-      vx  = lyapunov_symm(A,b1*b1');
-      vv = diag(aa*vx*aa'+b2*b2');
-      for i=1:M_.exo_nbr
-	vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)');
-	gamma_y{nar+2}(:,i) = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'))./vv;
-      end
-    end
-  else
-    lambda = options_.hp_filter;
-    ngrid = options_.hp_ngrid;
-    freqs = 0 : ((2*pi)/ngrid) : (2*pi*(1 - .5/ngrid)); 
-    tpos  = exp( sqrt(-1)*freqs);
-    tneg  =  exp(-sqrt(-1)*freqs);
-    hp1 = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2);
-    
-    mathp_col = [];
-    IA = eye(size(A,1));
-    IE = eye(M_.exo_nbr);
-    for ig = 1:ngrid
-      f_omega  =(1/(2*pi))*( [inv(IA-A*tneg(ig))*ghu1;IE]...
-			     *M_.Sigma_e*[ghu1'*inv(IA-A'*tpos(ig)) ...
-		    IE]); % state variables
-      g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % selected variables
-      f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
-      mathp_col = [mathp_col ; (f_hp(:))'];    % store as matrix row
-                                               % for ifft
-    end;
-
-    % covariance of filtered series
-    imathp_col = real(ifft(mathp_col))*(2*pi);
-
-    gamma_y{1} = reshape(imathp_col(1,:),nvar,nvar);
-    
-    % autocorrelations
-    if nar > 0
-      sy = sqrt(diag(gamma_y{1}));
-      sy = sy *sy';
-      for i=1:nar
-	gamma_y{i+1} = reshape(imathp_col(i+1,:),nvar,nvar)./sy;
-      end
-    end
-    
-    %variance decomposition
-    if M_.exo_nbr > 1 
-      gamma_y{nar+2} = zeros(nvar,M_.exo_nbr);
-      SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord)=M_.Sigma_e+1e-14*eye(M_.exo_nbr);
-      cs = chol(SS)';
-      SS = cs*cs';
-      b1(:,M_.exo_names_orig_ord) = ghu1;
-      b2(:,M_.exo_names_orig_ord) = ghu(iky,:);
-      mathp_col = [];
-      IA = eye(size(A,1));
-      IE = eye(M_.exo_nbr);
-      for ig = 1:ngrid
-	f_omega  =(1/(2*pi))*( [inv(IA-A*tneg(ig))*b1;IE]...
-			       *SS*[b1'*inv(IA-A'*tpos(ig)) ...
-		    IE]); % state variables
-	g_omega = [aa*tneg(ig) b2]*f_omega*[aa'*tpos(ig); b2']; % selected variables
-	f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
-	mathp_col = [mathp_col ; (f_hp(:))'];    % store as matrix row
-						 % for ifft
-      end;
-
-      imathp_col = real(ifft(mathp_col))*(2*pi);
-      vv = diag(reshape(imathp_col(1,:),nvar,nvar));
-      for i=1:M_.exo_nbr
-	mathp_col = [];
-	SSi = cs(:,i)*cs(:,i)';
-	for ig = 1:ngrid
-	  f_omega  =(1/(2*pi))*( [inv(IA-A*tneg(ig))*b1;IE]...
-				 *SSi*[b1'*inv(IA-A'*tpos(ig)) ...
-		    IE]); % state variables
-	  g_omega = [aa*tneg(ig) b2]*f_omega*[aa'*tpos(ig); b2']; % selected variables
-	  f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
-	  mathp_col = [mathp_col ; (f_hp(:))'];    % store as matrix row
-						   % for ifft
-	end;
-
-	imathp_col = real(ifft(mathp_col))*(2*pi);
-	gamma_y{nar+2}(:,i) = abs(diag(reshape(imathp_col(1,:),nvar,nvar)))./vv;
-      end
-    end
-  end
-  if sscanf(version('-release'),'%d') < 13
-    warning on
-  else
-    eval('warning on MATLAB:dividebyzero')
-  end
-  
-  % 10/18/2002 MJ
-  % 10/30/2002 added autocorrelations, gamma_y is now a cell array
-  % 01/20/2003 MJ added variance decomposition
-  % 02/18/2003 MJ added HP filtering (Thanks to Jean Chateau for the code)
-  % 04/28/2003 MJ changed handling of options
-  % 05/19/2003 MJ don't assume lags are in increasing order in building A
-  % 05/21/2003 MJ added global i_exo_var_,
-  %               variance decomposition: test M_.exo_nbr > 1
-  % 05/29/2003 MJ removed global i_exo_var_
-  % 06/10/2003 MJ test release for warning syntax
\ No newline at end of file
+% Copyright (C) 2001 Michel Juillard
+%
+% computes the theoretical auto-covariances, Gamma_y, for an AR(p) process 
+% with coefficients dr.ghx and dr.ghu and shock variances Sigma_e_
+% for a subset of variables ivar (indices in lgy_)
+% Theoretical HP filtering is available as an option
+
+function [Gamma_y,ivar]=th_autocovariances(dr,ivar)
+  global M_ options_
+
+  exo_names_orig_ord  = M_.exo_names_orig_ord;
+  if sscanf(version('-release'),'%d') < 13
+    warning off
+  else
+    eval('warning off MATLAB:dividebyzero')
+  end
+  nar = options_.ar;
+  Gamma_y = cell(nar+1,1);
+  if isempty(ivar)
+    ivar = [1:M_.endo_nbr]';
+  end
+  nvar = size(ivar,1);
+  
+  ghx = dr.ghx;
+  ghu = dr.ghu;
+  npred = dr.npred;
+  nstatic = dr.nstatic;
+  kstate = dr.kstate;
+  order = dr.order_var;
+  iv(order) = [1:length(order)];
+  nx = size(ghx,2);
+  
+  ikx = [nstatic+1:nstatic+npred];
+  
+  A = zeros(nx,nx);
+  k0 = kstate(find(kstate(:,2) <= M_.maximum_lag+1),:);
+  i0 = find(k0(:,2) == M_.maximum_lag+1);
+  i00 = i0;
+  n0 = length(i0);
+  A(i0,:) = ghx(ikx,:);
+  AS = ghx(:,i0);
+  ghu1 = zeros(nx,M_.exo_nbr);
+  ghu1(i0,:) = ghu(ikx,:);
+  for i=M_.maximum_lag:-1:2
+    i1 = find(k0(:,2) == i);
+    n1 = size(i1,1);
+    j1 = zeros(n1,1);
+    j2 = j1;
+    for k1 = 1:n1
+      j1(k1) = find(k0(i00,1)==k0(i1(k1),1));
+      j2(k1) = find(k0(i0,1)==k0(i1(k1),1));
+    end
+    AS(:,j1) = AS(:,j1)+ghx(:,i1);
+    i0 = i1;
+  end
+  b = ghu1*M_.Sigma_e*ghu1';
+
+
+  [A,B] = kalman_transition_matrix(dr);
+  % index of predetermined variables in A
+  i_pred = [nstatic+(1:npred) M_.endo_nbr+1:length(A)];
+  A = A(i_pred,i_pred);
+
+  if options_.order == 2
+    [vx,ns_var] =  lyapunov_symm(A,b);
+    i_ivar = find(~ismember(ivar,dr.order_var(ns_var+nstatic)));
+    ivar = ivar(i_ivar);
+    iky = iv(ivar);  
+    aa = ghx(iky,:);
+    bb = ghu(iky,:);
+    Ex = (dr.ghs2(ikx)+dr.ghxx(ikx,:)*vx(:)+dr.ghuu(ikx,:)*M_.Sigma_e(:))/2;
+    Ex = (eye(n0)-AS(ikx,:))\Ex;
+    Gamma_y{nar+3} = AS(iky,:)*Ex+(dr.ghs2(iky)+dr.ghxx(iky,:)*vx(:)+dr.ghuu(iky,:)*M_.Sigma_e(:))/2;
+  end
+  if options_.hp_filter == 0
+    if options_.order < 2
+      [vx, ns_var] =  lyapunov_symm(A,b);
+      i_ivar = find(~ismember(ivar,dr.order_var(ns_var+nstatic)));
+      ivar = ivar(i_ivar);
+      iky = iv(ivar);  
+      aa = ghx(iky,:);
+      bb = ghu(iky,:);
+    end
+    Gamma_y{1} = aa*vx*aa'+ bb*M_.Sigma_e*bb';
+    k = find(abs(Gamma_y{1}) < 1e-12);
+    Gamma_y{1}(k) = 0;
+    
+    % autocorrelations
+    if nar > 0
+      vxy = (A*vx*aa'+ghu1*M_.Sigma_e*bb');
+
+      sy = sqrt(diag(Gamma_y{1}));
+      sy = sy *sy';
+      Gamma_y{2} = aa*vxy./sy;
+      
+      for i=2:nar
+	vxy = A*vxy;
+	Gamma_y{i+1} = aa*vxy./sy;
+      end
+    end
+    
+    % variance decomposition
+    if M_.exo_nbr > 1
+      Gamma_y{nar+2} = zeros(length(ivar),M_.exo_nbr);
+      SS(exo_names_orig_ord,exo_names_orig_ord)=M_.Sigma_e+1e-14*eye(M_.exo_nbr);
+      cs = chol(SS)';
+      b1(:,exo_names_orig_ord) = ghu1;
+      b1 = b1*cs;
+      b2(:,exo_names_orig_ord) = ghu(iky,:);
+      b2 = b2*cs;
+      vx  = lyapunov_symm(A,b1*b1');
+      vv = diag(aa*vx*aa'+b2*b2');
+      for i=1:M_.exo_nbr
+	vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)');
+	Gamma_y{nar+2}(:,i) = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'))./vv;
+      end
+    end
+  else
+    if options_.order < 2
+      iky = iv(ivar);  
+      aa = ghx(iky,:);
+      bb = ghu(iky,:);
+    end
+    lambda = options_.hp_filter;
+    ngrid = options_.hp_ngrid;
+    freqs = 0 : ((2*pi)/ngrid) : (2*pi*(1 - .5/ngrid)); 
+    tpos  = exp( sqrt(-1)*freqs);
+    tneg  =  exp(-sqrt(-1)*freqs);
+    hp1 = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2);
+    
+    mathp_col = [];
+    IA = eye(size(A,1));
+    IE = eye(M_.exo_nbr);
+    for ig = 1:ngrid
+      f_omega  =(1/(2*pi))*( [inv(IA-A*tneg(ig))*ghu1;IE]...
+			     *M_.Sigma_e*[ghu1'*inv(IA-A'*tpos(ig)) ...
+		    IE]); % state variables
+      g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % selected variables
+      f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
+      mathp_col = [mathp_col ; (f_hp(:))'];    % store as matrix row
+                                               % for ifft
+    end;
+
+    % covariance of filtered series
+    imathp_col = real(ifft(mathp_col))*(2*pi);
+
+    Gamma_y{1} = reshape(imathp_col(1,:),nvar,nvar);
+    
+    % autocorrelations
+    if nar > 0
+      sy = sqrt(diag(Gamma_y{1}));
+      sy = sy *sy';
+      for i=1:nar
+	Gamma_y{i+1} = reshape(imathp_col(i+1,:),nvar,nvar)./sy;
+      end
+    end
+    
+    %variance decomposition
+    if M_.exo_nbr > 1 
+      Gamma_y{nar+2} = zeros(nvar,M_.exo_nbr);
+      SS(exo_names_orig_ord,exo_names_orig_ord)=M_.Sigma_e+1e-14*eye(M_.exo_nbr);
+      cs = chol(SS)';
+      SS = cs*cs';
+      b1(:,exo_names_orig_ord) = ghu1;
+      b2(:,exo_names_orig_ord) = ghu(iky,:);
+      mathp_col = [];
+      IA = eye(size(A,1));
+      IE = eye(M_.exo_nbr);
+      for ig = 1:ngrid
+	f_omega  =(1/(2*pi))*( [inv(IA-A*tneg(ig))*b1;IE]...
+			       *SS*[b1'*inv(IA-A'*tpos(ig)) ...
+		    IE]); % state variables
+	g_omega = [aa*tneg(ig) b2]*f_omega*[aa'*tpos(ig); b2']; % selected variables
+	f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
+	mathp_col = [mathp_col ; (f_hp(:))'];    % store as matrix row
+						 % for ifft
+      end;
+
+      imathp_col = real(ifft(mathp_col))*(2*pi);
+      vv = diag(reshape(imathp_col(1,:),nvar,nvar));
+      for i=1:M_.exo_nbr
+	mathp_col = [];
+	SSi = cs(:,i)*cs(:,i)';
+	for ig = 1:ngrid
+	  f_omega  =(1/(2*pi))*( [inv(IA-A*tneg(ig))*b1;IE]...
+				 *SSi*[b1'*inv(IA-A'*tpos(ig)) ...
+		    IE]); % state variables
+	  g_omega = [aa*tneg(ig) b2]*f_omega*[aa'*tpos(ig); b2']; % selected variables
+	  f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
+	  mathp_col = [mathp_col ; (f_hp(:))'];    % store as matrix row
+						   % for ifft
+	end;
+
+	imathp_col = real(ifft(mathp_col))*(2*pi);
+	Gamma_y{nar+2}(:,i) = abs(diag(reshape(imathp_col(1,:),nvar,nvar)))./vv;
+      end
+    end
+  end
+  if sscanf(version('-release'),'%d') < 13
+    warning on
+  else
+    eval('warning on MATLAB:dividebyzero')
+  end
+  
\ No newline at end of file