# ------------------------------------------------------ # Program: TimeSeries3.R # Author: Steven M. Boker # Date: Tue Nov 7 10:04:59 EST 2006 # # Timeseries, State Space and Vector Plots Part 3 # # ------------------------------------------------------ # ------------------------------------------------ # Load required libraries. library(lattice) # ------------------------------------------------------ # Create a cyclic time series and # synchronous cross-regression relation tX <- sin(seq(0,20, by=.1)) * .7 tY <- 2 + .5 * tX[1:201] + rnorm(201,mean=0,sd=.2) tX <- tX[1:201] + rnorm(201,mean=0,sd=.2) tData <- data.frame(x=tX,y=tY) summary(tData) pdf("TimeSeries3XY1.pdf", height=5 ,width=6) plot(c(1,length(tX)), c(-1,3), type="n" , xlab="Sample", ylab="Value") lines(c(1:length(tX)), tX, type="l", col="blue") lines(c(1:length(tX)), tY, type="l", col="red") dev.off() # ------------------------------------------------------ # Create an cross-correlation function plot of tX and tY maxTau <- 100 tLen <- length(tX) tCCF1 <- rep(NA, maxTau+1) for (tau in (-maxTau:maxTau)) { if (tau >= 0) { t1 <- tX[1:(tLen-tau)] t2 <- tY[(1+tau):tLen] } else { t1 <- tY[1:(tLen+tau)] t2 <- tX[(1-tau):tLen] } tSelect <- !is.na(t1) & !is.na(t2) if (length(t1[tSelect]) < 5) next tCCF1[maxTau+tau+1] <- cor(t1[tSelect], t2[tSelect]) } pdf("TimeSeries3Cross1.pdf", height=5 ,width=6) plot(c(-maxTau,maxTau), c(-1,1), type="n" , xlab="Lag", ylab="Cross-correlation") lines(c(-maxTau:maxTau), tCCF1, type="l", lty=1, lwd=2, col=4) lines(c(-maxTau,maxTau), c(0,0), type="l", lty=2) lines(c(0,0), c(-1,1), type="l", lty=2) dev.off() # ------------------------------------------------------ # Create a cyclic time series and # time-delayed cross-regression relation tX <- sin(seq(0,20, by=.1)) * .7 tY <- 2 + .5 * tX[21:201] + rnorm(181,mean=0,sd=.2) tX <- tX[1:181] + rnorm(181,mean=0,sd=.2) tData <- data.frame(x=tX,y=tY) summary(tData) pdf("TimeSeries3XY2.pdf", height=5 ,width=6) plot(c(1,length(tX)), c(-1,3), type="n" , xlab="Sample", ylab="Value") lines(c(1:length(tX)), tX, type="l", col="blue") lines(c(1:length(tX)), tY, type="l", col="red") dev.off() # ------------------------------------------------------ # Create an cross-correlation function plot of tX and tY maxTau <- 100 tLen <- length(tX) tCCF1 <- rep(NA, maxTau+1) for (tau in (-maxTau:maxTau)) { if (tau >= 0) { t1 <- tX[1:(tLen-tau)] t2 <- tY[(1+tau):tLen] } else { t1 <- tY[1:(tLen+tau)] t2 <- tX[(1-tau):tLen] } tSelect <- !is.na(t1) & !is.na(t2) if (length(t1[tSelect]) < 5) next tCCF1[maxTau+tau+1] <- cor(t1[tSelect], t2[tSelect]) } pdf("TimeSeries3Cross2.pdf", height=5 ,width=6) plot(c(-maxTau,maxTau), c(-1,1), type="n" , xlab="Lag", ylab="Cross-correlation") lines(c(-maxTau:maxTau), tCCF1, type="l", lty=1, lwd=2, col=4) lines(c(-maxTau,maxTau), c(0,0), type="l", lty=2) lines(c(0,0), c(-1,1), type="l", lty=2) dev.off() # ------------------------------------------------------ # Write the timeseries matrix to a textfile . write(tY, "CrossSine2.dat", ncolumns=1) # ------------------------------------------------------ # Create a randomized sequence independent variable, a regression # relation and a "Problem" condition. rX <- rep(c(1,2,3,4),50)[order(runif(200,0,1))] rY <- 1 + .8*rX + rnorm(200,mean=0,sd=.5) tSelect1 <- c(F,rX[1:199]==4) tSelect2 <- c(rX[1:199]==4,F) rY[tSelect1] <- rY[tSelect2] rDataFrame <- data.frame(rX, rY) summary(rDataFrame) # ------------------------------------------------------ # Calculate a regression summary(lm(rY~rX, data=rDataFrame)) # ------------------------------------------------------ # Scatterplot the data pdf("TimeSeries3Plot1.pdf", height=5 ,width=6) plot(c(0,5), c(0,6), type="n" , xlab="X", ylab="Y") lines(rX, rY, type="p", col="blue") abline(lm(rY~rX, data=rDataFrame)) dev.off() # ------------------------------------------------------ # Create a timelagged dataframe tau <- 1 tLen <- length(rX) rDataLagged <- data.frame(x1=rX[1:(tLen-tau)], x2=rX[(1+tau):tLen], y1=rY[1:(tLen-tau)], y2=rY[(1+tau):tLen]) # ------------------------------------------------------ # Calculate first-order autocorrelations and crosscorrelations round(cor(rDataLagged),2) # --------------------------------------------------------------------- # Pairs Scatterplot pdf("TimeSeries3Pairs1.pdf", height=5 ,width=5) pairs(rDataLagged, labels=c("X(t)", "X(t+1)", "Y(t)", "Y(t+1)"), cex=.5, col="blue") dev.off() # --------------------------------------------------------------------- # Pairs Scatterplot selecting the problem area. tSelect <- rDataLagged\$y1 > 3.5 & rDataLagged\$y2 > 3.5 pdf("TimeSeries3Pairs1b.pdf", height=5 ,width=5) pairs(as.matrix(rDataLagged)[tSelect,], labels=c("X(t)", "X(t+1)", "Y(t)", "Y(t+1)"), cex=.5, col="blue") dev.off() # --------------------------------------------------------------------- # Pairs Scatterplot selecting X1==4. tSelect <- rDataLagged\$x1==4 pdf("TimeSeries3Pairs1c.pdf", height=5 ,width=5) pairs(as.matrix(rDataLagged)[tSelect,], labels=c("X(t)", "X(t+1)", "Y(t)", "Y(t+1)"), cex=.5, col="blue") dev.off() # ------------------------------------------------------ # Write the timeseries matrix to a textfile . tMatrix <- as.matrix(rDataLagged) write(round(t(tMatrix),4)[3,], "rDataLagged.dat", ncolumns=dim(tMatrix)[2]) # ------------------------------------------------------ # Create a timelagged dataframe (tau = 2) tau <- 2 tLen <- length(rX) rDataLagged2 <- data.frame(x1=rX[1:(tLen-tau)], x2=rX[(1+tau):tLen], y1=rY[1:(tLen-tau)], y2=rY[(1+tau):tLen]) # ------------------------------------------------------ # Calculate first-order autocorrelations and crosscorrelations round(cor(rDataLagged2),2) # --------------------------------------------------------------------- # Pairs Scatterplot pdf("TimeSeries3Pairs2.pdf", height=5 ,width=5) pairs(rDataLagged2, labels=c("X(t)", "X(t+2)", "Y(t)", "Y(t+2)"), cex=.5, col="blue") dev.off() # ------------------------------------------------------ # Create a cyclic time series and # synchronous cross-regression relation tX <- sin(seq(0,20, by=.1)) tY <- 5 + .5 * tX[1:201] + rnorm(201,mean=0,sd=.05) tX <- tX[1:201] + rnorm(201,mean=0,sd=.05) tData <- data.frame(x=tX,y=tY) summary(tData) # ------------------------------------------------------ # Create a timelagged dataframe (3-D Embedding) tau <- 15 tLen <- length(tX) rDataLagged3D <- data.frame(x1=tX[1:(tLen-(2*tau))], x2=tX[(1+tau):(tLen-tau)], x3=tX[(1+(2*tau)):tLen], y1=tY[1:(tLen-(2*tau))], y2=tY[(1+tau):(tLen-tau)], y3=tY[(1+(2*tau)):tLen]) for(i in seq(0,360,by=10)) { tfilename <- paste("TimeSeries3ThreeD-", i, ".pdf", sep="") pdf(tfilename, height=5 ,width=6) print(cloud(x1 ~ x2 * x3, data = rDataLagged3D, screen = list(z = i, x = -70, y = 3))) dev.off() } # ------------------------------------------------------ # Quit the program