%%Written by Erin Boykin, 2011
%————————————————————————————————-
%Inputs: X – time series for analysis (rows are time point, columns are
% channels)
% p – max order of AR model fit (perform an AIC or BIC analysis prior
% in order to make sure max is set sufficiently high)
% freq – vector of discrete frequencies at which gc is calculated
%
%Outputs: gc – gc values in form frequency x channel i x channel j x
% channel k, where causality is from channel i to channel j
% conditioned on channel k
%————————————————————————–
%The ARfit package for computing AR models is available here : http://www.gps.caltech.edu/~tapio/arfit/
% it is required in order for this function to run
%
%Also note that this method uses the Partition Matrix approach. For more
%details refer to Chen et al. J Neurosci Methods 150, 228-37.
function [gc] = conditionalGC(X,p,freq)
Nc = size(X,2);
gc=zeros(length(freq),Nc,Nc,Nc);%Initialize all causality to zero
%First generate all arfit matrices, this will improve performance later!
all2D = cell(Nc,Nc);
all2Sig = cell(Nc,Nc);
for xx = 1:Nc
for zz = xx+1:Nc
if xx~=zz
[w,D,all2Sig{xx,zz}]=arfit(X(:,[xx zz]),0,p,’fpe’,'zero’);
all2D{xx,zz} = getBf(2,D,freq);
all2D{zz,xx} = all2D{xx,zz}([2 1],[2 1],:);
all2Sig{zz,xx} = all2Sig{xx,zz}([2 1],[2 1]);
end
end
end
all3D = cell(Nc,Nc,Nc);
all3Sig = cell(Nc,Nc,Nc);
for xx = 1:Nc
for zz = xx+1:Nc
for yy = zz+1:Nc
[w,D,all3Sig{xx,yy,zz}]=arfit(X(:,[xx yy zz]),0,p,’fpe’,'zero’);
all3D{xx,yy,zz} = getBf(3,D,freq);
all3D{xx,zz,yy} = all3D{xx,yy,zz}([1 3 2],[1 3 2],:);
all3D{yy,xx,zz} = all3D{xx,yy,zz}([2 1 3],[2 1 3],:);
all3D{yy,zz,xx} = all3D{xx,yy,zz}([2 3 1],[2 3 1],:);
all3D{zz,xx,yy} = all3D{xx,yy,zz}([3 1 2],[3 1 2],:);
all3D{zz,yy,xx} = all3D{xx,yy,zz}([3 2 1],[3 2 1],:);
all3Sig{xx,zz,yy} = all3Sig{xx,yy,zz}([1 3 2],[1 3 2]);
all3Sig{yy,xx,zz} = all3Sig{xx,yy,zz}([2 1 3],[2 1 3]);
all3Sig{yy,zz,xx} = all3Sig{xx,yy,zz}([2 3 1],[2 3 1]);
all3Sig{zz,xx,yy} = all3Sig{xx,yy,zz}([3 1 2],[3 1 2]);
all3Sig{zz,yy,xx} = all3Sig{xx,yy,zz}([3 2 1],[3 2 1]);
end
end
end
for xx = 1:Nc
for zz = 1:Nc
if xx~=zz
%%Consider case where X is driven, Z is condtioned on
%Do Geweke’s 2D transform
all2D{xx,zz}(2,1,:) = all2D{xx,zz}(2,1,:) + (all2Sig{xx,zz}(1,2)/all2Sig{xx,zz}(1,1))*all2D{xx,zz}(1,1,:);
all2D{xx,zz}(1,1,:) = all2D{xx,zz}(1,1,:) + (all2Sig{xx,zz}(1,2)/all2Sig{xx,zz}(1,1))*all2D{xx,zz}(1,2,:);
all2Sig{xx,zz}(2,2) = all2Sig{xx,zz}(2,2) – all2Sig{xx,zz}(1,2)^2/all2Sig{xx,zz}(1,1);
for yy = 1:Nc
if yy~=xx && yy~=zz
%Do Geweke’s 3D transform
H = getGewekeTransform3(all3D{xx,yy,zz},all3Sig{xx,yy,zz});
[SigHat,G] = partitionMethod(H,all3Sig{xx,yy,zz});
Gpad = [G(1,1,:),zeros(1,1,length(freq)),G(1,2,:);zeros(1,1,length(freq)),ones(1,1,length(freq)),zeros(1,1,length(freq));G(2,1,:),zeros(1,1,length(freq)),G(2,2,:)];
for ff = 1:length(freq)
Q(:,:,ff) = (Gpad(:,:,ff)^(-1))*H(:,:,ff);
gc(ff,yy,xx,zz) = log(abs(SigHat(1,1,ff))/abs(Q(1,1,ff)*all3Sig{xx,yy,zz}(1,1)*Q(1,1,ff)’));
end
end
end
end
end
end
function [Bf] = getBf(m,B,freq);
popt = size(B,2) / m;
A0 = eye(m);
Bf = zeros(m,m,length(freq));
for g=1:length(freq)
for ch=1:m
Bj = B(:,[ch:m:m*popt]);
for j = 1:popt
Bf(:,ch,g) = Bf(:,ch,g) + Bj(:,j).*exp(1i*j*freq(g));
end
end
Bf(:,:,g) = A0 – Bf(:,:,g);
end
function [SigHat,G] = partitionMethod(H,Sig)
HHat = zeros(2,1,size(H,3));
SigHat = zeros(2,2,size(H,3));
P = zeros(2,2,size(H,3));
G = zeros(2,2,size(H,3));
for ff=1:size(H,3)
HHat(:,:,ff) = H([1 3],[1 3],ff)^-1*H([1 3],[2],ff);
SigHat(:,:,ff) = Sig([1 3],[1 3]) + HHat(:,:,ff)*[Sig(1,2),Sig(3,2)] + [Sig(1,2);Sig(3,2)]*[HHat(1,1,ff)',HHat(2,1,ff)'] + Sig(2,2)*HHat(:,:,ff)*[HHat(1,1,ff)',HHat(2,1,ff)'];
P(:,:,ff) = eye(2);
P(2,1,ff) = -1*(SigHat(1,2,ff)/SigHat(1,1,ff));
G(:,:,ff) = H([1 3],[1 3],ff)*P(:,:,ff)^-1;
end
function [HHat] = getGewekeTransform3(H,Sig)
P1 = eye(3);
P1(2,1) = -Sig(2,1)/Sig(1,1);
P1(3,1) = -Sig(3,1)/Sig(1,1);
P2 = eye(3);
P2(3,2) = -(Sig(3,2)-Sig(3,1)*Sig(1,2)/Sig(1,1))/(Sig(2,2)-Sig(2,1)*Sig(1,2)/Sig(1,1));
P = P2*P1;
HHat = zeros(size(H));
for ff = 1:size(H,3)
HHat(:,:,ff) = H(:,:,ff)*P^-1;
end
