# --------------------------------------------------------------------- # 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 # # --------------------------------------------------------------------- # ------------------------------------------------------ # Create an error-free autoregressive time series and # synchronous cross-regression relation maxSamples <- 201 tX1 <- rep(NA, maxSamples) tX1[1] <- 90 for (i in 2:maxSamples) { tX1[i] <- .975 * tX1[i-1] } tY1 <- 10 + .5 * tX1 + rnorm(maxSamples, mean=0, sd=5) tX1 <- 10 + tX1 + rnorm(maxSamples, mean=0, sd=5) tX1[tX1 <= 0] <- 0.001 # ------------------------------------------------------ # Create a timeseries dataframe of tX1 and tY1 tTime <- seq(0,1000, by=5)/60 tData1 <- data.frame(x=tX1, y=tY1, time=tTime) # --------------------------------------------------------------------- # Plot the two timeseries pdf("Decay1x.pdf", height=5, width=6) plot(tData1$time, tData1$x, xlab="Time", ylab="X", type='l', pch=1) dev.off() pdf("Decay1y.pdf", height=5, width=6) plot(tData1$time, tData1$y, xlab="Time", ylab="Y", type='l', pch=1) dev.off() # ----------------------------- # Create data vectors to hold autocorrelations and cross correlations totalLags <- 50 irXX <- rep(0, totalLags) irYY <- rep(0, totalLags) irYX <- rep(0, totalLags) irXY <- rep(0, totalLags) # ----------------------------- # Calculate the autocorrelations and cross correlations tLength <- length(tData1$x) for (tau in 1:totalLags) { irXX[tau] <- cor(tData1$x[1:(tLength-(tau-1))], tData1$x[tau:tLength]) irYY[tau] <- cor(tData1$y[1:(tLength-(tau-1))], tData1$y[tau:tLength]) irYX[tau] <- cor(tData1$y[1:(tLength-(tau-1))], tData1$x[tau:tLength]) irXY[tau] <- cor(tData1$x[1:(tLength-(tau-1))], tData1$y[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("Decay1xAutoCorr.pdf", height=5, width=6) plot(c(-totalLags,totalLags), c(-1,1), xlab="Tau", ylab="Correlation", main="Autocorrelation Decay 1 X", type="n") lines(seq(-totalLags,totalLags),irXX2[1:((totalLags*2) + 1)], type="l") dev.off() pdf("Decay1yAutoCorr.pdf", height=5, width=6) plot(c(-totalLags,totalLags), c(-1,1), xlab="Tau", ylab="Correlation", main="Autocorrelation Decay 1 Y", type="n") lines(seq(-totalLags,totalLags),irYY2[1:((totalLags*2) + 1)], type="l") dev.off() pdf("Decay1xyCrossCorr.pdf", height=5, width=6) plot(c(-totalLags,totalLags), c(-1,1), xlab="Tau", ylab="Correlation", main="Cross Correlation Decay 1 X and Y", type="n") lines(seq(-totalLags,totalLags),irYX2[1:((totalLags*2) + 1)], type="l") dev.off() # --------------------------------------------------------------------- # Plot state space plots of X and Y theTau <- 1 tLen <- length(tData1$x) pdf("Decay1xStateSpaceTau1.pdf", height=5, width=5) plot(tData1$x[1:(tLen-theTau)], tData1$x[(1+theTau):tLen], xlab="x(t)", ylab="X(t+1)", type='l', pch=1) dev.off() pdf("Decay1yStateSpaceTau1.pdf", height=5, width=5) plot(tData1$y[1:(tLen-theTau)], tData1$y[(1+theTau):tLen], xlab="y(t)", ylab="y(t+1)", type='l', pch=1) dev.off() # --------------------------------------------------------------------- # --------------------------------------------------------------------- # --------------------------------------------------------------------- # --------------------------------------------------------------------- # Create an error dependent autoregressive time series and # timelagged cross-regression relation tX2 <- rep(NA, maxSamples+20) tX2[1] <- 1 for (i in 2:(maxSamples+20)) { tX2[i] <- .9 * tX2[i-1] + rnorm(1, mean=0, sd=2) } tY2 <- 5 + .5 * tX2[21:(maxSamples+20)] + rnorm(maxSamples, mean=0, sd=3) tX2 <- tX2[1:maxSamples] + rnorm(maxSamples, mean=0, sd=1) # ------------------------------------------------------ # Create a timeseries dataframe of tX2 and tY2 tTime <- seq(0,1000, by=5)/60 tData2 <- data.frame(x=tX2, y=tY2, time=tTime) # --------------------------------------------------------------------- # Plot the two timeseries pdf("Decay2x.pdf", height=5, width=6) plot(tData1$time, tData2$x, xlab="Time", ylab="X", type='l', pch=1) dev.off() pdf("Decay2y.pdf", height=5, width=6) plot(tData1$time, tData2$y, xlab="Time", ylab="Y", type='l', pch=1) dev.off() # ----------------------------- # Create data vectors to hold autocorrelations and cross correlations totalLags <- 50 irXX <- rep(0, totalLags) irYY <- rep(0, totalLags) irYX <- rep(0, totalLags) irXY <- rep(0, totalLags) # ----------------------------- # Calculate the autocorrelations and cross correlations tLength <- length(tData2$x) for (tau in 1:totalLags) { irXX[tau] <- cor(tData2$x[1:(tLength-(tau-1))], tData2$x[tau:tLength]) irYY[tau] <- cor(tData2$y[1:(tLength-(tau-1))], tData2$y[tau:tLength]) irYX[tau] <- cor(tData2$y[1:(tLength-(tau-1))], tData2$x[tau:tLength]) irXY[tau] <- cor(tData2$x[1:(tLength-(tau-1))], tData2$y[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("Decay2xAutoCorr.pdf", height=5, width=6) plot(c(-totalLags,totalLags), c(-1,1), xlab="Tau", ylab="Correlation", main="Autocorrelation Decay 2 X", type="n") lines(seq(-totalLags,totalLags),irXX2[1:((totalLags*2) + 1)], type="l") dev.off() pdf("Decay2yAutoCorr.pdf", height=5, width=6) plot(c(-totalLags,totalLags), c(-1,1), xlab="Tau", ylab="Correlation", main="Autocorrelation Decay 2 Y", type="n") lines(seq(-totalLags,totalLags),irYY2[1:((totalLags*2) + 1)], type="l") dev.off() pdf("Decay2xyCrossCorr.pdf", height=5, width=6) plot(c(-totalLags,totalLags), c(-1,1), xlab="Tau", ylab="Correlation", main="Cross Correlation Decay 2 X and Y", type="n") lines(seq(-totalLags,totalLags),irYX2[1:((totalLags*2) + 1)], type="l") dev.off() # --------------------------------------------------------------------- # Plot state space plots of X and Y theTau <- 1 tLen <- length(tData2$x) pdf("Decay2xStateSpaceTau1.pdf", height=5, width=5) plot(tData2$x[1:(tLen-theTau)], tData2$x[(1+theTau):tLen], xlab="x(t)", ylab="X(t+1)", type='l', pch=1) dev.off() pdf("Decay2yStateSpaceTau1.pdf", height=5, width=5) plot(tData2$y[1:(tLen-theTau)], tData2$y[(1+theTau):tLen], xlab="y(t)", ylab="y(t+1)", type='l', pch=1) dev.off() # --------------------------------------------------------------------- # --------------------------------------------------------------------- # --------------------------------------------------------------------- # --------------------------------------------------------------------- # Simulate two composite sine waves maxSamples <- 1000 # ----------------------------- # Pick two frequencies at random theFreq1 <- runif(1, min=.1, max=.5) theFreq2 <- runif(1, min=.05, max=.1) # ----------------------------- # 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 <- 20 tSine2 <- 0.5 * tSine1[(1+phaseOffset):maxSamples] + rnorm((maxSamples-phaseOffset), mean=0, sd=1) # --------------------------------------------------------------------- # Plot the two composite sine waves pdf("tSine1.pdf", height=5, width=6) plot(c(0,500), c(-2,2), xlab="Samples", ylab="Displacement", main="Composite Sine Wave 1", type="n") lines(seq(1,500), tSine1[1:500], type="l") dev.off() pdf("tSine2.pdf", height=5, width=6) plot(c(0,500), c(-2,2), xlab="Samples", ylab="Displacement", main="Composite Sine Wave 2", type="n") lines(seq(1,500), tSine2[1:500], type="l") dev.off() # ----------------------------- # Create data vectors to hold autocorrelations and cross correlations totalLags <- 100 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,totalLags),irXX2[1:201], 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,totalLags),irYY2[1:201], 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,totalLags),irYX2[1:201], type="l") dev.off() # --------------------------------------------------------------------- # End the program q()