# How to obtain cross-correlation function from cross-spectrum?

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?

## 1 Answer

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) #  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.220189 # variance inferred from the mapping Re(b) #  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")) 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.

• 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). – user68144 Mar 7 '15 at 18:14
• Probably the number of observations in the irregularly spaced series, but I'm not familiar with the periodogram for unevenly sampled time series. – javlacalle Mar 7 '15 at 22:38