[R] Spectral Coherence

Joseph Park jpark.us at att.net
Tue Jul 12 13:09:14 CEST 2011

   Thanks you! I should have realized that without explicitly engaging
   some form of averaging (which raises a windowing question) that the
   coh is always 1.
   On 7/12/2011 4:48 AM, Rolf Turner wrote:

     On 12/07/11 09:04, Joseph Park wrote:

         I would like to estimate a spectral coherence between
         two timeseries. The stats : spectrum() returns a coh matrix
         which estimates coherence (squared).
         A basic test which from which i expect near-zero coherence:
         x = rnorm(500)
         y = rnorm(500)
         xts = ts(x, frequency = 10)
         yts = ts(y, frequency = 10)
         gxy = spectrum( cbind( xts, yts ) )
         plot( gxy $ freq, gxy $ coh )
         yields a white spectrum of 1. Clearly i'm not using
         this correctly... or i mis-interpret the coh as a cross-spectral
         density estimate of coherence |Gxy|^2/(Gxx Gyy)
         Thanks in advance!

     By default spectrum() calls spec.pgram() with spans=NULL.  The result is
     it calculates the coherence between x and y as
             |I_{xy}(omega)|^2 / (I_{xx}(omega) * I_{yy}(omega))
     where I_{xy}() is the cross periodogram and I_{xx}() and I_{yy}() are the
     respective periodograms.
     This quantity will indeed be identically 1 --- see equation 10.1 in
     Bloomfield, second ed.,
     page 203.
     It would be nice if the help on spectrum mentioned this.
     According to Bloomfield what is needed is not the periodograms but rather
     estimated spectra,
     smoothed  versions of the periodograms, s_{xy}() etc.  Equation 10.4 in
     Bloomfield, second
     ed.,  page  206  indicates  that  the spectral estimates satisfy an
         |s_{xy}(omega)|^2 <= s_{xx}(omega) * s_{yy}(omega)
     whence the coherence is always between 0 and 1.
     To resolve the problem you need to specify the spans argument in the call
     to spectrum,
         gxy <- spectrum(cbind(xts,yts),spans=c(5,7))
             Rolf Turner

More information about the R-help mailing list