function [tab,tab1]=vardec(A,S,nz,C)
%VARDEC   Compute variance decompositions
%         [tab,tab1]=vardec(A,S,nz,C)
%

m  = length(nz);
n  = length(S);
Vx = dlyap(A,S);
V  = diag(Vx);
for i=1:m;
  k=int2str(nz(i));
  eval(['i',k,'=zeros(n);']) 
  eval(['i',k,'(',k,',',k,')=1;']) 
end;
s = chol(S)';
tab = [];
for j=1:m;
  k=int2str(nz(j));
  eval(['vx',k,'=dlyap(A,s*i',k,'*s'');'])
  eval(['v',k,'=diag(vx',k,');'])
  eval(['tab=[tab,v',k,'./V];'])
end;

if nargin > 3;
  W = diag(C*Vx*C');
  tab1 = [];
  for j=1:m;
    k=int2str(nz(j));
    eval(['v',k,'=diag(C*vx',k,'*C'');'])
    eval(['tab1=[tab1,v',k,'./W];'])
  end;
end;

