docs: def granger(vec1,vec2,order=10,rate=200,maxfreq=0):
      """
      GRANGER
      
      Provide a simple way of calculating the key quantities.
      
      Usage:
          F,pp,cohe,Fx2y,Fy2x,Fxy=granger(vec1,vec2,order,rate,maxfreq)
      where:
          F is a 1xN vector of frequencies
          pp is a 2xN array of power spectra
          cohe is the coherence between vec1 and vec2
          Fx2y is the causality from vec1->vec2
          Fy2x is the causality from vec2->vec1
          Fxy is non-directional causality (cohe-Fx2y-Fy2x)
          
          vec1 is a time series of length N
          vec2 is another time series of length N
          rate is the sampling rate, in Hz
          maxfreq is the maximum frequency to be returned, in Hz
      
      Version: 2011jul18
      """
      from .bsmart import timefreq, pwcausalr
      from scipy import array, size
      
      if maxfreq==0: F=timefreq(vec1,rate) # Define the frequency points
      else: F=array(list(range(0,maxfreq+1))) # Or just pick them
      npts=size(F,0)
      
      data=array([vec1,vec2])
      F,pp,cohe,Fx2y,Fy2x,Fxy=pwcausalr(data,1,npts,order,rate,maxfreq)
      return F,pp[0,:],cohe[0,:],Fx2y[0,:],Fy2x[0,:],Fxy[0,:]