function [gp,se]=varci(A,X,p,nz,C,type);
%VARCI   Compute confidence intervals for variance decompositions of VAR's.
%        [v,V,int]=varci(Theta,W,func,type) assumes that Theta are estimates
%        underlying A,Sigma in 
%       
%                 x[t+1] = A x[t] + eps[t+1],  eps ~ N(0,Sigma).
%
%        that are normally distributed, N(Theta,W). 
%        The mapping from Theta to A,Sigma is given in a .m file with
%        name "func".  The input func is a string. The type of method
%        used for obtaining the intervals depends on type: if type=1,
%        a normal approximation is used; if type=2, a bootstrap 
%        method is used.
%        [v,V,int,tab]=varci(Theta,func,type) produces a table, tab, with
%        the results. 
%
[T,n]=size(X);
XX=ones(T-p,1);
for i=1:p;
  XX=[XX,X(p-i+1:T-i,:)];
end;
p
y=X(p+1:T,:);
size(XX)
size(A')
size(y)
res=y-XX*A';
S=cov(res);
W=kron(S,inv(XX'*XX));

Theta=A(:);
[mm,nn]=size(A);
n=length(Theta);
m=length(nz);
del=diag(max(abs(Theta)*1e-5,1e-8));
for i=1:m;
  k=int2str(nz(i));
  eval(['i',k,'=zeros(S);']) 
  eval(['i',k,'(',k,',',k,')=1;']) 
end;
G=[];
a1=zeros(A);
a2=zeros(A);
for i=1:n;
  if rem(i,10)==0; disp(i); end;
  a1(:)=Theta+del(:,i);
  S=cov(y-XX*a1')+eye(S)*1e-13;
  V = diag(C*edlyap(a1(:,2:nn),S)*C');
  s = chol(S)';
  gp= [];
  for j=1:m;
    k=int2str(nz(j));
    eval(['v',k,'=diag(C*edlyap(a1(:,2:nn),s*i',k,'*s'')*C'');'])
    eval(['gp=[gp,v',k,'./V];'])
  end;
  gp=gp(:);
  a2(:)=Theta-del(:,i);
  S=cov(y-XX*a2')+eye(S)*1e-13;
  V = diag(C*edlyap(a2(:,2:nn),S)*C');
  s = chol(S)';
  gm= [];
  for j=1:m;
    k=int2str(nz(j));
    eval(['v',k,'=diag(C*edlyap(a2(:,2:nn),s*i',k,'*s'')*C'');'])
    eval(['gm=[gm,v',k,'./V];'])
  end;
  gm=gm(:);
  G = [G,(gp-gm)/(2*del(i,i))];
end;
save stuff0

S=cov(y-XX*A')+eye(S)*1e-13;
V = diag(C*edlyap(A(:,2:nn),S)*C');
s = chol(S)';
gp= [];
for j=1:m;
  k=int2str(nz(j));
  eval(['v',k,'=diag(C*edlyap(A(:,2:nn),s*i',k,'*s'')*C'');'])
  eval(['gp=[gp,v',k,'./V];'])
end;

gp

size(G)
Vg = G*W*G'
se = zeros(length(V),m);
se(:) = diag(sqrt(Vg));      
se

