function [v,V,int,tab]=varci(Theta,W,func,nz,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. 
%

n=length(Theta);
m=length(nz);
del=diag(max(abs(Theta)*1e-1,1e-8));
eval(['[A,S]=',func,'(Theta);'])
for i=1:m;
  k=int2str(nz(i));
  eval(['i',k,'=zeros(S);']) 
  eval(['i',k,'(',k,',',k,')=1;']) 
end;
G=[];
for i=1:n;
  eval(['[A,S]=',func,'(Theta+del(:,i));'])
  V = diag(edlyap(A,S));
  s = chol(S)';
  gp= [];
  for j=1:m;
    k=int2str(nz(j));
    eval(['v',k,'=diag(edlyap(A,s*i',k,'*s''));'])
    eval(['gp=[gp,v',k,'./V];'])
  end;
  gp=gp(:);
  eval(['[A,S]=',func,'(Theta-del(:,i));'])
  V = diag(edlyap(A,S));
  s = chol(S)';
  gm= [];
  for j=1:m;
    k=int2str(nz(j));
    eval(['v',k,'=diag(edlyap(A,s*i',k,'*s''));'])
    eval(['gm=[gm,v',k,'./V];'])
  end;
  gm=gm(:);
  G = [G,(gp-gm)/(2*del(i,i))];
end;
G
Vg = G*W*G'
se = zeros(length(V),m);
se(:) = diag(sqrt(Vg));      

eval(['[A,S]=',func,'(Theta);'])
V = diag(edlyap(A,S));
s = chol(S)';
gp= [];
for j=1:m;
  k=int2str(nz(j));
  eval(['v',k,'=diag(edlyap(A,s*i',k,'*s''));'])
  eval(['gp=[gp,v',k,'./V];'])
end;

gp
se



