function [tabtot,tabind]=evardec(A,B,C,J,K,Q)
%VARDEC calculates variance decompositions of the j-step ahead prediction 
%  errors, v(j)=E[z(t+j)-E[z(t+j) | x(t)] ][z(t+j)-E[z(t+j) | x(t)] ]', 
%  j in J, for the system:
%
%      x(t+1) = A x(t) + B w(t)
%        z(t) = C x(t) 
%
%      E [w(t) w(t)'] = Q
%
%  For this system, v(j) = C*B*Q*B'*C' + C*A*B*Q*B'*A'*C' + ...
%                                      + C*A^(j-1)*B*Q*B'*A^(j-1)'*C',  j>0
%                       =v(j-1) + C*A^(j-1)*B*Q*B'*A^(j-1)'*C'
%  where v(0)=0.
%  To calculate the j-step ahead variance due to the kth component of w(t),
%  k in K, we compute vk(j)=E[z(t+j)-E[z(t+j) | w(k,t),w(k,t-1),...]]*
%  [z(t+j)-E[z(t+j) | w(k,t),w(k,t-1),...]]'.  To do this, we use a Cholesky
%  decomposition: Q=R'*R and let D=B*R'.  So, B*w(t)=D*u(t), where Eu(t)*
%  u(t)'=I.  Then v(j)=C*D*D'*C'+ ... +C*A^(j-1)*D*D'*A^(j-1)'*C' and
%  vk(j)=C*D*i(k)*D'*C'+...+C*A^(j-1)*D*i(k)*D'*A^(j-1)'*C', where i(k)
%  is a matrix of zeros and element (k,k) equal 1.
%
%  INPUTS:  A, B, C, J, K, Q  (or just A, B, C, J, K,  if Q=I)
%    Notes:
%      If J=-N, for some N, then I compute v(j), vk(j) for j=1,..N.
%      If K=0, I don't compute vk(j)
%      If K=-1, I compute vk(j) for all possible k.
%      Otherwise, I use the j's k's in J, K as specified
%
%  OUTPUTS: 
%    (i)  tabtot = a table with a row equal to diag(v(j)) for each j in J.
%    (ii) tabind = a table giving error by error results: diag(vk(j)), j in J
%                  and k in K.
%    Note: for each table, a column with elements of J appears before columns of
%          prediction error variances.
%
ra=max(size(A));
[rc,col]=size(C);
rc1=rc+1;
if nargin==5; Q=eye(ra); end;
R=chol(Q);
D=B*R';
if min(J)<0; J=[0:-min(J)]'; end; 
if K==-1; K=[1:ra]'; end;
lJ=length(J);
if K~=0; lK=length(K); else; lK=0; end;
N=max(J);
sort(J);
sort(K);
cc=0;
vj=zeros(rc);
vjk=zeros(rc,lK*rc);
for j=1:N;
  AI=A^(j-1);
  vj =vj  + C*AI*D*D'*AI'*C';
  for k=1:lK;
    ik=zeros(ra);
    ik(K(k),K(k))=1;
    vjk(:,rc*(k-1)+1:rc*k)=vjk(:,rc*(k-1)+1:rc*k)+ C*AI*D*ik*D'*AI'*C';
  end;
  if min(abs(j-J))==0;
    cc=cc+1;
    tabtot(cc,:)=[j,diag(vj)'];
    for k=1:lK;
      tabind(cc,rc1*(k-1)+1:rc1*k)=[j,diag(vjk(:,rc*(k-1)+1:rc*k))'];
    end;
  end;
end;
if K==0; tabind=0; end;
