Hello MIke,

Thank you for taking the time to answer this question.

I noticed in your book and the online code accompanying chapter 28, you do not discuss nor show how to calculate Granger Prediction in the frequency domain at the trial level.

Looking at your online code, Figure 28.5, you first reshape a 3d matrix of ( electrode by time by trial ) into a 2d matrix (electrode by time*trial). After this, you run armorf() function to get your autoregression error values.

I have slightly modified your code below which is essentially unchanged, however, I loop through trials rather than combining into a single long trial. Your book is so comprehensive and I am curious that there is a particular reason that this was left out. As far as you can tell, is the below script methodologically appropriate?

```
tempdata = eegdata(2, times, trialn);
```

for triali = 1:size(tempdata,3)

[Ax,Ex] = armorf(tempdata(1,:, triali), 1, window_length, order_points);

[Ay,Ey] = armorf(tempdata(2,:, triali), 1, window_length, order_points);

[Axy,E] = armorf(tempdata , 1, window_length, order_points);

% fit AR models (model estimation from bsmart toolbox)

% this is the original code, to be run on reshaped data

%[Ax,Ex] = armorf(tempdata(1,:), trialn, window_length, order_points);

%[Ay,Ey] = armorf(tempdata(2,:), trialn , window_length, order_points);

%[Axy,E] = armorf(tempdata , trialn, window_length, order_points);

```
% code below is adapted from bsmart toolbox function pwcausal.m
% corrected covariance
eyx = E(2,2) - E(1,2)^2/E(1,1);
exy = E(1,1) - E(2,1)^2/E(2,2);
N = size(E,1);
for fi=1:length(frequencies)
% transfer matrix (note the similarity to Fourier transform)
H = eye(N);
for m = 1:order_points
H = H + Axy(:,(m-1)*N+1:m*N)*exp(-1i*m*2*pi*frequencies(fi)/srate);
end
Hi = inv(H);
S = H\E*Hi'/srate;
% granger prediction per frequency
```

%mini is the sliding time window index at which I store these results.

```
tf_granger(1,fi,mini, triali) = log( abs(S(2,2))/abs(S(2,2)-(Hi(2,1)*exy*conj(Hi(2,1)))/srate) );
tf_granger(2,fi,mini, triali) = log( abs(S(1,1))/abs(S(1,1)-(Hi(1,2)*eyx*conj(Hi(1,2)))/srate) );
end
```

end

Thank you so much!

Hiro