# --------------------------------------------------------------------- # Program: ACF-CCF-1.R # Author: Steven M. Boker # Date: Thu Feb 22 10:10:31 EST 2007 # # This simulates two composites and then calculates autocorrelation # and cross correlation and then plots the results. # # # --------------------------------------------------------------------- # Revision History # Steven M. Boker -- Thu Feb 22 10:11:23 EST 2007 # Created ACF-CCF-1.R from Autocorrelation1.R # # --------------------------------------------------------------------- # --------------------------------------------------------------------- # --------------------------------------------------------------------- # --------------------------------------------------------------------- # --------------------------------------------------------------------- # Simulate a composite sine waves maxSamples <- 300 # ----------------------------- # Pick two frequencies at random #theFreq1 <- runif(1, min=.1, max=.5) #theFreq2 <- runif(1, min=.05, max=.1) theFreq1 <- .3 # ----------------------------- # The first data vector is a composite of the two frequencies tSine1 <- sin(theFreq1*(1:maxSamples)) # + sin(theFreq2*(1:maxSamples)) # ----------------------------- # The second data vector is dependent on the first time series phaseOffset <- 0 tSine2 <- 0.5 * tSine1[(1+phaseOffset):maxSamples] + rnorm((maxSamples-phaseOffset), mean=0, sd=.75) # --------------------------------------------------------------------- # Plot the two composite sine waves pdf("tSine1.pdf", height=5, width=6) plot(c(0,100), c(-2,2), xlab="Samples", ylab="Displacement", main="Composite Sine Wave 1", type="n") lines(seq(1,200), tSine1[1:200], type="l") dev.off() pdf("tSine2.pdf", height=5, width=6) plot(c(0,100), c(-2,2), xlab="Samples", ylab="Displacement", main="Composite Sine Wave 2", type="n") lines(seq(1,200), tSine2[1:200], type="l") dev.off() # ----------------------------- # Create data vectors to hold autocorrelations and cross correlations totalLags <- 30 irXX <- rep(0, totalLags) irYY <- rep(0, totalLags) irYX <- rep(0, totalLags) irXY <- rep(0, totalLags) # ----------------------------- # Calculate the autocorrelations and cross correlations tLength <- length(tSine2) for (tau in 1:totalLags) { irXX[tau] <- cor(tSine1[1:(tLength-(tau-1))],tSine1[tau:tLength]) irYY[tau] <- cor(tSine2[1:(tLength-(tau-1))],tSine2[tau:tLength]) irYX[tau] <- cor(tSine2[1:(tLength-(tau-1))],tSine1[tau:tLength]) irXY[tau] <- cor(tSine1[1:(tLength-(tau-1))],tSine2[tau:tLength]) } # ----------------------------- # Paste together the leads and lags revSeq <- seq(totalLags, 2, -1) irXX2 <- c(irXX[revSeq],irXX) irYX2 <- c(irXY[revSeq],irYX) irYY2 <- c(irYY[revSeq],irYY) # --------------------------------------------------------------------- # Plot the autocorrelation and cross correlation between the data vectors pdf("tSine1AutoCorr.pdf", height=5, width=6) plot(c(-totalLags,totalLags), c(-1,1), xlab="Tau", ylab="Correlation", main="Autocorrelation Sine Wave 1", type="n") lines(seq(-totalLags+1,totalLags+1),irXX2[1:((totalLags*2) + 1)], type="l") lines(c(-totalLags+1,totalLags+1), c(0,0), type="l") dev.off() pdf("tSine2AutoCorr.pdf", height=5, width=6) plot(c(-totalLags,totalLags), c(-1,1), xlab="Tau", ylab="Correlation", main="Autocorrelation Sine Wave 2", type="n") lines(seq(-totalLags+1,totalLags+1),irYY2[1:((totalLags*2) + 1)], type="l") lines(c(-totalLags+1,totalLags+1), c(0,0), type="l") dev.off() pdf("tSine12CrossCorr.pdf", height=5, width=6) plot(c(-totalLags,totalLags), c(-1,1), xlab="Tau", ylab="Correlation", main="Cross Correlation Sine Waves 1 and 2", type="n") lines(seq(-totalLags+1,totalLags+1),irYX2[1:((totalLags*2) + 1)], type="l") lines(c(-totalLags+1,totalLags+1), c(0,0), type="l") dev.off() # --------------------------------------------------------------------- # Plot state space plots of X and Y tLen <- length(tSine2) for (theTau in 1:8) { pdf(paste("tSine1StateSpaceTau", theTau, ".pdf", sep=""), height=5, width=5) plot(tSine1[1:(tLen-theTau)], tSine1[(1+theTau):tLen], xlab="x(t)", ylab=paste("x(t+", theTau, ")", sep=""), type='p', pch=1) dev.off() pdf(paste("tSine2StateSpaceTau", theTau, ".pdf", sep=""), height=5, width=5) plot(tSine2[1:(tLen-theTau)], tSine2[(1+theTau):tLen], xlab="y(t)", ylab=paste("y(t+", theTau, ")", sep=""), type='p', pch=1) dev.off() } # ----------------------------------------------------------------------- # Create 100 time-shuffled surrogates of the Lorenz data # and calculate their mean ACF maxSurrogates <- 1000 surrogateACFMatrix <- matrix(NA, maxSurrogates, totalLags) tLength <- length(tSine2) for (i in 1:maxSurrogates) { permutedSeries <- tSine2[order(runif(tLength))] for (tau in 1:totalLags) { surrogateACFMatrix[i, tau] <- cor(permutedSeries[1:(tLength-(tau-1))],permutedSeries[tau:tLength]) } } surrogateACF985 <- rep(NA, totalLags) surrogateACF025 <- rep(NA, totalLags) for (tau in 1:totalLags) { surrogateACF985[tau] <- sort(surrogateACFMatrix[, tau])[985] surrogateACF025[tau] <- sort(surrogateACFMatrix[, tau])[025] } surrogateACF985B <- c(surrogateACF985[seq(totalLags, 2, -1)], surrogateACF985) surrogateACF025B <- c(surrogateACF025[seq(totalLags, 2, -1)], surrogateACF025) # ----------------------------------------------------------------------- # Plot ACF of tSine2 with critical values pdf("tSine2ACFwithCritValues.pdf", height=5, width=6) plot(c(-totalLags,totalLags), c(-1,1), xlab="Tau", ylab="Correlation", main="Autocorrelation Sine Wave 2", type="n") lines(seq(-totalLags+1,totalLags+1),irYY2[1:((totalLags*2) + 1)], type="l") lines(seq(-totalLags+1,totalLags+1),surrogateACF985B[1:((totalLags*2) + 1)], type="l", lty=2) lines(seq(-totalLags+1,totalLags+1),surrogateACF025B[1:((totalLags*2) + 1)], type="l", lty=2) lines(c(-totalLags+1,totalLags+1), c(0,0), type="l") dev.off() # --------------------------------------------------------------------- # End the program #q()