pygc/
├── parametric.py
├── non_parametric.py
├── granger.py
├── tools.py
├── pySpec.py

import pygc.parametric
import pygc.granger
import pygc.pySpec

X = # Data which you want to estimate GC, dimension must be [Nvariables, Nobservations, Ntrials]
Nvars  = X.shape[0]  # For now I recommend Nvars to be equal two.
N      = X.shape[1]
Trials = X.shape[2]
Fs = # Sample frequency of the data

# Computing the frequency axis
f = pygc.pySpec.compute_freq(N, Fs)

m   = 2 # Model order you wish to fit
AR  = np.zeros([m,Nvars,Nvars])  # Store the auto-regressive coefficients
SIG = np.zeros([Nvars,Nvars])    # Store the noise covariance matrix
for T in range(Trials):
  aux1, aux2 = pygc.parametric.YuleWalker(X[:,:,T], m, maxlags=100)
	AR  += aux1/Trials
	SIG += aux2/Trials

# The below works for Nvars = 2, for more variables you should do this in a pairwise fashion
# Compute transfer matrix (H) and spectral matrix (S) both with size [Nvars, Nvars,len(f)]
H, S = pygc.parametric.compute_transfer_function(AR, SIG, f, Fs)
# Return the three components of GC
Ix2y, Iy2x, Ixy  = pygc.granger.granger_causality(S, H, SIG)

import pygc.non_parametric
import pygc.granger
import pygc.pySpec

X = # Data which you want to estimate GC, dimension must be [Nvariables, Nobservations, Ntrials]
Nvars  = X.shape[0]  # For now I recommend Nvars to be equal two.
N      = X.shape[1]
Trials = X.shape[2]
Fs = # Sample frequency of the data

# Computing the frequency axis
f = pygc.pySpec.compute_freq(N, Fs)

# Declare spectral matrix
S = np.zeros([Nvars,Nvars,N//2+1]) + 1j*np.zeros([Nvars,Nvars,N//2+1])
# Compute the avarage spectral matrix for each Trial
# Here I use FFT to compute the spectrum you can use other methods you like
for i in range(Trials):
	S[0,0] += pygc.pySpec.cxy(X=X[0,:,T], Y=[], f=f, Fs=Fs) / Trials
	S[0,1] += pygc.pySpec.cxy(X=X[0,:,T], Y=Z[1,i,:], f=f, Fs=Fs) / Trials
	S[1,0] += pygc.pySpec.cxy(X=X[1,:,T], Y=Z[0,i,:], f=f, Fs=Fs) / Trials
	S[1,1] += pygc.pySpec.cxy(X=X[1,:,T], Y=[], f=f, Fs=Fs) / Trials

# The below works for Nvars = 2, for more variables you should do this in a pairwise fashion
# Uses Wilson factorization to estimate the transfer function and the noise covariance matrix
Snew, H, SIG = pygc.non_parametric.wilson_factorization(S, f, Fs, Niterations=30)
# Return the three components of GC
Ix2y, Iy2x, Ixy  = pygc.granger.granger_causality(S, H, SIG)

# Now Nvars >= 3
S = np.zeros([Nvars, Nvars, N//2 + 1]) * (1 + 1j)

# Compute the spectral matrix
for trial in range(Trials):
	for i in range(Nvars):
		for j in range(nvars):
			S[i,j] += pygc.pySpec.cxy(X=Y[i,:,trial], Y=Y[j,:,trial], f=f, Fs=Fs) / Trials

# Estimate the conditional GC in time domain
F   = pygc.granger.conditional_granger_causality(S, f, Fs, Niterations = 30, verbose=False)
# Estimate the conditional GC in frequency domain
cGC = pygc.granger.conditional_spec_granger_causality(S, f, Fs, Niterations=30, tol=1e-12, verbose=True)