1
$\begingroup$

I have implemented a piece of code based on the Lomb-Scargle approach for determining the cross-spectrum of two time series. My cross spectrum contains complex numbers and I have used the basic fft function in R with the option inverse=true to apply a inverse Fast Fourier Transform to obtain the cross-correlation but i get a complex valued cross correlation function which is not correct. What other methods should I use?

$\endgroup$
2
$\begingroup$

I will give you an answer for the mapping between the autocovariances and the spectrum of a time series. Your question is about the cross-correlation and cross-spectrum, but the ideas below may be helpful.

The mapping between the autocovariances or order $\tau$, $\gamma(\tau)$ and the spectral density, $f(\omega)$, can be expressed as follows:

$$ \gamma(\tau) = \int_{-\pi}^\pi e^{i\omega\tau} f(\omega) \, d\omega \,. $$

Given the spectral density, the autocovariances can be obtained in R as follows:

# generate a sample series
set.seed(123)
x <- arima.sim(n=200, model=list(ar=c(0.6,-0.3,0.4)))
# sample autocovariances (to be revored from the periodogram)
gamma <- acf(x, type="cov", plot = FALSE)$acf[,,1]
# spectral density estimate (smoothed periodogram)
sde <- spectrum(x, spans = c(3,3), plot=FALSE)
# apply the mapping in the equation above for the frequencies at which 
# the spectral density was estimated
w <- sde$freq # frequencies (omega)
b <- rep(NA, length(gamma))
for (tau in seq.int(0, length(gamma)-1))
{
  tmp <- exp(1i * tau * w * 2 * pi) * sde$spec
  b[tau+1] <- sum(tmp + Conj(tmp))/length(x)
}
head(cbind(gamma, Re(b)))
# [1,] 1.21408770 1.27297299
# [2,] 0.47336471 0.50364341
# [3,] 0.03438318 0.05592229
# [4,] 0.30775323 0.36804157
# [5,] 0.28104297 0.32918129
# [6,] 0.09878623 0.11652643
tail(cbind(gamma, Re(b)))
# [19,] -0.07139406 -0.087051153
# [20,] -0.14771390 -0.158636206
# [21,] -0.19801861 -0.194581525
# [22,] -0.07825159 -0.075446230
# [23,]  0.02632276 -0.004064849
# [24,] -0.10278970 -0.105784687
# Imaginary terms cancel out
Im(b)
# [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# sample variance
var(x)
# [1] 1.220189
# variance inferred from the mapping
Re(b)[1]
# [1] 1.272973

Sample autocovariances and the values inferred from the mapping are close to each other:

plot(seq.int(0, length(gamma)-1), gamma, xlab = "lag order")
lines(seq.int(0, length(gamma)-1), Re(b))
legend("topright", bty = "n", lty = c(0, 1), pch = c(1, NA), legend =
  c("sample autocovariances", "autocovariances mapped from the spectral density"))

Mapping spectral density to autocovariances

The following points will also apply for the cross-spectrum and cross-covariances:

  • More accurate values are expected when a smoothed version of the periodogram (as obtained by spectrum) rather than the raw periogram Mod(fft(x))/(2*pi*length(x)).

  • Depending on the range and scale of frequencies used to compute the spectrum, your definition of $\omega$ (w in the code) may change.

  • Observe that in the code above the integral is summed over the range $(0, \pi]$. The conjugate values are then added up so that the complete range $(-\pi, \pi]$ is covered. The imaginary part of the complex numbers will therefore cancel out.

  • There may be a way to use the inverse Fourier transform that you mention, but for a test version it may be safer to obtain the products $e^{i\omega\tau} f(\omega)$ in an explicit loop.

$\endgroup$
  • $\begingroup$ I initially start with two time series of sample size equal to 100, then I sample them to have two irregularly sampled time series of mean sampling intervals 1 and 4. In this case, what would be the normalization factor at the last line( the one where you divide by the length of x). $\endgroup$ – user68144 Mar 7 '15 at 18:14
  • $\begingroup$ Probably the number of observations in the irregularly spaced series, but I'm not familiar with the periodogram for unevenly sampled time series. $\endgroup$ – javlacalle Mar 7 '15 at 22:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.