On 12/07/11 09:04, Joseph Park wrote:
Greetings,
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
that
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 *inequality*
|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,
e.g.
gxy <- spectrum(cbind(xts,yts),spans=c(5,7))
HTH
cheers,
Rolf Turner
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
